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ABSTRACT 


The military’s dependence on fossil fuels for electric power production in isolated settings 
is both logistically and monetarily expensive. Currently, the Department of Defense is ac¬ 
tively seeking alternative methods to produce electricity, thus decreasing dependence on 
fossil fuels and increasing combat power. We believe piezoelectric generators have the 
ability to contribute to military applications of alternative electrical power generation in 
isolated and austere conditions. In this paper, we use three and six variable mathematical 
models to analyze piezoelectric generator power generation capabilities. Using m k factorial 
sampling, nearly orthogonal and balanced Latin hypercube (NOBLH) design, and NOBLH 
iterative methods, we find optimal solutions to maximize piezoelectric generator power 
output. We further analyze our optimal results using robustness analysis techniques to de¬ 
termine the sensitivity of our models to variable precision. With our results, we provide 
analysts and engineers the optimal designs involving material parameters in the piezoelec¬ 
tric generator, as well as the generator’s environment, in order to maximize electric output. 
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CHAPTER 1: 
Introduction 


1.1 Motivation 

"Unleash us from the tether of fuel." 

-Lieutenant General Jamis Mattis, 

Future Fuels, Naval Research 
Advisory Committee Report 
April 2006 

During the past decade of conflict, service members on bases across Iraq and Afghanistan 
could hear the constant hum of diesel generators. These thirsty machines provide the base’s 
electric lifeblood enabling lighting, refrigeration, life-support heating and cooling, battery 
charging, and equipment operation. Years of conflict proved our dependence on fossil fuels 
comes with a cost. A constant consumption of diesel fuel requires a hefty logistics tail 
to provide fossil fuels to every base across the area of operations. During Operation Iraqi 
Freedom and Operation Enduring Freedom an Army study found that, for every 24 convoys 
sent out, one service member or civilian engaged in fuel transport was killed. Additionally, 
due to transportation cost of fossil fuels to isolated bases, the price of fuel escalates from 
approximately one dollar per gallon to, in some cases, 400 dollars per gallon [1]. 

The source of electricity generation service members depend on decreases combat capabil¬ 
ity. Military forces are taken away from combat operations and nation building to guard 
and operate convoys. Dollars spent on expensive fossil fuels are better allocated to support 
service members and operations. Additionally, expeditionary forces do not have the ability 
to transport heavy fuels. Their dependence on fossil fuels for power generation inhibits 
their ability to conduct operations. 

This paper attempts to improve upon an alternative method of producing electricity while 
in austere conditions without the necessity of fossils fuels and their accompanying logistics 
tail. 
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Figure 1.1: Soldier guarding fuel convoy after attack [2]. 


1.2 Piezoelectricity Definition and History 


The name “piezo” is derived from the Greek piezo or piezein, meaning “to squeeze or 
press” [3]. Piezoelectricity is pressure-driven electricity. Certain materials generate elec¬ 
tricity when under pressure (mechanical stress). These same materials change shape when 
electricity is applied to the material. A common use of the piezoelectricity phenomenon is 
demonstrated in quartz watches. A battery in the watch feeds electricity to a piezoelectric 
quartz crystal. The electricity induces a precise vibration in the quartz. The mechanisms in 
the watch then convert the quartz vibrations into a means of keeping accurate time. 
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Figure 1.2: Piezoelectric quartz crystal resonator used for time keeping [4], 


In 1880, two French physicists, Pierre and Jacques Curie, found that in certain materials, 
such as zinc blende, sodium chlorate, boracite, tourmaline, cane sugar, Rochelle salt, and 
quartz, mechanical stresses induce macroscopic polarization and hence the production of 
electric surface charges. The following year, Gabriel Lippmann predicted the converse 
effect from fundamental thermodynamic principles, that is, an imposed voltage produces 
mechanical deformations or strains of the material [5]. In 1882, the Curie brothers con¬ 
firmed the existence of the converse effect based on experimental observations [3]. 

Pierre Curie used the piezoelectric effect to measure charges emitted by radium, but piezo¬ 
electricity was not used for practical applications for several decades. During World War I, 
Paul Langevin developed piezoelectric crystal quartz transducer depth sounding devices to 
locate submerged vessels, primarily German submarines [6]. In the 1920s, developments 
in crystal resonators for the stabilization of oscillators launched the field of frequency of 
control. The introduction of quartz control drastically changed the way humans keep time. 
Today, piezoelectric applications include smart materials for vibration control, aerospace 
and astronautical applications of flexible surfaces and structures, sensors for robotic appli¬ 
cations, and novel applications for vibration reduction in sports equipment (tennis racquets 



and snowboards) [5]. Recently, scientists have focused research on using piezoelectric 
materials for energy generation [7]—[14]. 


1.3 How Piezoelectricity Works 

The separation of positive and negative electrical charges in many molecules results in 
what is known as a dipole moment. In piezoelectric crystals, when a mechanical stress is 
applied (the crystal being compressed, twisted, or pulled) the molecular dipole moments 
re-orient themselves. This re-orientation causes a variation in surface charge density and 
thus produces voltage. The effect is illustrated in Figure 1.3 for forces normal to the mate¬ 
rial. Conversely, when an electric field is applied across a piezoelectric medium, there is a 
slight change in the shape of the dipoles causing a very small, but significant change in the 
material dimensions. 

Pulling normal force Equilibrium Compressional normal force 

AAA 



V V V 

Figure 1.3: Illustration of the piezoelectric concept [3]. 
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1.4 Applications of Piezoelectric Power Generation 

Using the converse piezoelectric effect, it is possible to generate electricity by deforming 
piezoelectric materials. However, piezoelectric generation is far from the most efficient 
method of power generation. Presently, piezoelectric generators can only produce small 
amounts of electricity, but with the increase of efficiency of powered devices, piezoelectric 
power generation has become a possible supply of electricity for certain devices. 

Devices well suited for piezoelectric power supply include the following: 

• devices with low power requirements 

• devices without power supplies in close proximity 

• devices without easy access to a power supply (power lines are inconvenient or im¬ 
possible) 

• devices without access to other alternative power sources (solar, wind, geothermal) 

• devices with access to vibrations or other mechanical manipulation 

The most obvious device that meets the above criteria is a sensor. Sensors are often located 
far from power sources, concealed from light sources, and require power to collect and 
send information. It is not always possible or convenient to connect power to the sensors 
by wire. Examples of sensors that could receive power through piezoelectric generation 
include sensors in air ducts [15], tires, and inside machinery. 

Piezoelectric power generation has the potential to enhance military operations as well. 
Military units often operate in primitive/isolated environments where it is not practical to 
transport conventional fuel to generate electricity. Military units must rely on alternative 
methods of power generation or do without. In addition to solar and wind power genera¬ 
tion, service members could generate power with piezoelectric power generators fixed to 
personnel, equipment, or any other item in motion. 


1.5 Research Objectives 

This paper attempts to find the optimal designs involving material parameters in the piezo¬ 
electric generator, as well as the generator’s environment, in order to maximize electric 
output. 
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1.6 Thesis Organization 

This chapter introduced the motivation and relevant background of harvesting piezoelectric 
energy. In Chapter 2, we review a one-dimensional power generation model and study how 
to optimize power generation using variables in the model. From the results of Chapter 2 
we derive a three-variable power generation model. In Chapter 3, we derive a six-variable 
power generation model using a single mode model. In Chapter 4, we describe the various 
methods we use in order to optimize our three-variable and six-variable power generation 
models. In Chapter 5, we discuss the results of our optimization methods and the robust¬ 
ness of our results. In Chapter 6, we conclude the paper and discuss further studies to be 
completed based on the information and ideas generated in this paper. Finally, we include 
an appendix to display MATLAB code used for calculations within the paper. 
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CHAPTER 2: 

One-Dimensional Power Generation Model: 
Derivation of the Three Variable Model 


2.1 1-D Model Design: Interpretation Scheme 

We can use a basic 1-D model to analyze the power generated from a vibration energy 
harvester to understand conversion. This model is limited to harvesters where the electri¬ 
cal damping term is linear and proportional to velocity. Nevertheless, this simple model 
is useful in understanding the feasibility of the device and input parameters on power ex¬ 
tracted. The electrical energy is extracted from the mechanical system, which is excited 
by a mechanical input. This extraction is not necessarily linear, or proportional to velocity, 
however, it is a dissipative process and can generally be viewed as electrical damping. 

Following [16], we consider a generator (shown in Figure 2.1) which consists of a seismic 
mass m on a spring k. When the generator is vibrated, the mass moves out of phase with the 
generator housing, so that there is a net movement between the mass and the housing. This 
relative displacement is sinusoidal in amplitude, and can drive a suitable transducer to gen¬ 
erate electrical energy. The transducer is depicted as a dashpot, d , because the conversion 
of mechanical energy into electrical energy damps the mass. There are several transduction 
methods suitable for the generator. The choice of transducer makes little difference to the 
amount of electrical power that will be generated. Three possible transduction mechanisms 
are (1) piezoelectric: using piezoelectric material to convert strain in the spring into elec¬ 
tricity; (2) electromagnetic: a magnet attached to the mass induces a voltage in a coil as it 
moves; and (3) electrostatic: an electric arrangement with a permanent charge embedded 
in the mass induces a voltage on the plates of a capacitor as it moves. 
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Figure 2.1: A schematic diagram of the generator [17]. 


In the simplified model, we assume that the mass of the vibration source is much bigger 
than the mass of the seismic mass in the generator, and the vibration source is an infinite 
source of power. This implies that the vibration source is unaffected by the movement of 
the generator. 

The differential equation of motion can be derived from the dynamic forces on the mass: 

mz ( t ) + dz ( t ) + kz ( t ) = -my" ( t ) (2.1) 

where the generator housing is vibrated with a displacement y(t), the relative motion of the 
mass with respect to the housing is z(t), m is the seismic mass, d is the damping constant, 
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and k is the spring constant. The solution to this second-order ordinary differential equation 
(ODE) with constant coefficients consists of two parts, the complementary function, which 
is the solution of the corresponding homogeneous equation, and the particular solution. The 
complementary function, in this case, is a damped free vibration which goes to zero as time 
increases. The particular solution is the one of interest here. For a sinusoidal excitational 
vibration, y(t) = Yocos(cot), the particular solution to the preceding ODE is a steady-state 
oscillation of the same frequency, co, as that of the excitation. We can assume the particular 
solution to be of the form 

z(t) = Acos(cot — 0) (2.2) 

where A is the amplitude of oscillation and (j) is the phase of the displacement with respect 
to the exciting force. 

Differentiating Equation (2.2) yields 


z' (t) — —Aft) sin (tut — (j)), 
z" (t) = —Aft) 2 cos (cot — (j)). 

Substituting the above equations into Equation (2.1) we have 


(2.3) 

(2.4) 


— Amcoi 2 cos (cot — (j)) — Adco sin (cot — (j)) +Akc os (cot — (j)) — mYoCO 2 cos (cot) . (2.5) 


Utilizing the trigonometric identities cos(cot — <j>) — cos(cot)cos((j)) + sin((Dt)sin(ty) and 
sin(cot — 0) = sin(cot)cos(<f>) — cos(cot)sin(<j>) we can rewrite Equation (2.5) as 


—Am oo 2 cos cot cos (f) — Am ft) 2 sin cursing — Ad 00 sin cot cos <f> + Ad co cos cot sin (j) 
+Ak cos cot cos (j> +Aksincotsin(t) = mY^co 2 cos (cot). 


( 2 . 6 ) 


Matching the coefficients for cos (cot) and sin (cot), respectively, we find 


(Ak — Amco 2 ) cos (j) +AJft)sin0 = mYoCO 2 , 

(2.7) 

(k — mco 2 ) sin (j) — dcocoscj) = 0. 

(2.8) 


9 



This leads to 


tan (j) — 


dco 


k — mco 


2 ’ 


[ (k — mco 2 ) cos (j) + dca sin 0] A — mYoCO 2 , 


[ (k — mco 2 ) + dco tan A = mYoCO 2 sec (j) = mYoCO 2 \J 1 + tan 2 ^. 

Solving for A, we obtain 


(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 


mY{)CD 2 \J 1 + tan 2 </> 
(k — mco 2 ) +dcotan(f) 


mY{)C0 2 \j\ + (^^i) 

(i k-mco 2 )+dco(- k 

mYoCO 2 

A = / = 

y (k — mco 2 ) 2 + (Jco) 2 


mYoCO 2 \J (k — mco 2 ) 2 + (Jft))^ 
(k — mco 2 ) 2 + (dco) 2 


( 2 . 12 ) 

(2.13) 


We now express the amplitude and phase in non-dimensional form. Dividing the numerator 
and denominator by fc, we get 


tan0 


dco 

k — mco 2 


dco 

k 


mco 2 5 
k 


(2.14) 


A = 



(2.15) 


These expressions can be further expressed in terms of the following quantities where CO n 
is the natural frequency of undamped oscillation (resonant angular frequency), d c is critical 
damping, and £ is the modal damping ratio: 


(On — 



(2.16) 


d c = 2nico n . 


dco 

k 


d d c CO 
d c k 


2 C’ 

1 


(2.17) 

(2.18) 
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(2.19) 


k k 
d c 2 m(O n 


The non-dimensional formulas for the amplitude and phase then become 


k/m 1 k (O n 

0 f~k 2 V m 2 


2C— 

tan(j) = -^ 


1 - -® 
V co„ 


( 2 . 20 ) 




( 2 . 21 ) 


The instantaneous power transfer in the mass is the product of the force on the mass and its 
velocity: 

P (t) = -my" (; t ) [/ (. t ) + z' (f)] . (2.22) 

When damping is present, due to the electrical transducer, there is a net transfer of mechan¬ 
ical power into electrical power. This net electrical power generated, P, is 

T T 

P = — Jp(t)dt — — j mco 2 Yocos (cot) [—YqCO sin (cot)—Acosin (cot — (j))\dt, (2.23) 


1 -T cos (2 cot) 


P— — J mAco Yo sin({)cos (cot)dt — — j mAco Yo sin $ -— 

0 0 


dt , (2.24) 


P = mYoCO 3 — sin (j). 

Substituting in Equation (2.21) for A, we obtain: 


(2.25) 


P = mY 0 2 CO 3 - 




2 a/1 T COt~(j) 


(2.26) 
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P = mY 0 2 (o 3 - 


1 - 5 ) +( 2 ^ 



(2.27) 


l-f 

\ (On J 

2C — 

b co„ 


P = 


mYo 2 co 3 


CO \ CO 

co„ / ^ ro~ 






1 


5) + ( 2 ^ 


i 


5) + ( 2 ^ 


2 ‘ 


(2.28) 


We now have an explicit expression for net electric power generated we can conduct analy¬ 
sis on. In the following sections we first view the output power as a function of normalized 
frequency and investigate how the output power behaves as we vary normalized frequency 
for various values of the modal damping ratio. After that, we treat the output power as 
a multi-variable function of both normalized frequency and modal damping ratio, seeking 
the optimal value of the output power. Then, we compare our model prediction with an 
experimental result. 


2.2 Power as a Function of Normalized Frequency 

First, we will treat net electric power generated, P , as a function of normalized frequency 
only. To do this, we will introduce C = m^Y^co 3 and x — to get the function 


P(x) 


C 


(1 —x 2 ) 2 + (2£.v) 


2' 


(2.29) 


In order to optimize the function we differentiate P(x) with respect to x: 


2 Cx 5 [x 4 — 4x 2 + 8£ 2 x 2 + 3] 
(. x 4 + 4£ 2 .v 2 — 2.v 2 + l) 2 


(2.30) 


By inspection, we can see if £ = 0, P'(x) does not exist when x = 1. In fact, P(l) becomes 
infinitely large. If £ ^ 0, the critical points of P(x) depend on the value of £. 
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2.2.1 Initial Analysis 

We will fix values for £ and find the value of x that maximizes P(x). Like Williams and 
Yates [17], we will analyze P(x) with £ values of 0.1,0.3,0.5, and 0.7. 

For £ = 0.1, 


2Cx 5 [x 4 -3.92x 2 + 3] 
(x 4 — 1.96.r 2 + l) 2 


(2.31) 


Solving P'{x) — 0 yields critical points x = ±1.6963, ±1.0211. Recall that x — > 0. For 

£ = 0.1, we plot (^(7 as a function of x in an interval containing the critical points in Figure 

2 . 2 . 



From Figure 2.2, we can see that the maximum power can be generated when x — j 4 - — 

iX)n 

1.0211 for £ = 0.1, and the maximum value of is 26.0421. 
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For C = 0.3, 


. 2C.r 5 \x 4 — 3.28.r 2 + 3] 

p w = - k - r 1 - 

(x 4 — l.64x 2 + 1) 

Solving P'(x) — 0 yields four complex roots. The critical point is x 
for £ = 0.3 in Figure 2.3. 


(2.32) 

0. Now we plot 



Figure 2.3 shows when £ = 0.3, is an increasing function of x. 
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For C = 0.5, 


P\x) 


2 Cx 5 [.x 4 — 2x 2 + 3] 
(x 4 — x 2 + 1)“ 


(2.33) 


Solving P'(x) — 0 yields four complex roots. The critical point is x = 0. Now we plot 
for £ = 0.5 in Figure 2.4. 



Figure 2.4 shows when £ = 0.5, ^4 is an increasing function of x. 
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For C = 0.7, 


. 2C.r 5 \.x 4 — 0.08x 2 + 3l 

p' ( x ) — _11 

(x 4 — 0.04x 2 + l) 2 

Solving P'(x) — 0 yields four complex roots. The critical point is x 
for C, — 0.7 in Figure 2.5 


(2.34) 

0. Now we plot 



Figure 2.5 shows when £ = 0.7, is an increasing function of x. 


16 





















2.2.2 Analysis as £ Approaches Zero 

One will notice only £ = 0.1 resulted in a critical point other than x — 0. Let us evaluate 
the critical points and corresponding optical when £ <0.1. For £ = 0.05, 


2Cx 5 [x 4 — 3.98.r 2 + 3] 
(x 4 -l.99x 2 + lf 


(2.35) 


Solving P'(x) = 0 yields critical points x — ±1.7233,3=1.0051. For £ = 0.05, we plot 
as a function of x in an interval containing the critical points in Figure 2.6. 



From Figure 2.6, we can see that the maximum power can be generated when x = -fjj- = 
1.0051 for £ = 0.05, and the maximum value of is 101.01. 
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It appears as £ approaches zero, the optimal value of jc approaches one and the resulting 

P(x) 

maximum power generation, -^ L , increases. See Figures 2.7 and 2.8 for further illustration 
of the how the optimal x value and maximum power generation change as £ decreases. 



£ 


Figure 2.7: Optimal x as a function of £. 


Figure 2.7 further depicts the optimal frequency ratio (^-) approaching one as £ approaches 
zero. This follows the principal of resonance, where a vibration is able to drive a sys¬ 
tem into larger oscillations when the vibration frequency matches the system’s natural fre¬ 
quency ((0 — co n ). 
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2500 



P(x) 

Figure 2.8 demonstrates that when we set a = (O n , -j 1 depends solely on how small we 
can set £. 

2.3 Power as a Function of Modal Damping Ratio and 
Normalized Frequency 

We will now treat net electric power generated, P, as a function of both modal damping 
ratio, C, and frequency ratio, x = jy. We can graphically solve for the optimal value of 
by plotting as the dependent variable and using £ and x—-^-as the independent 
variables. If we use the parameters 0.1 <£< 1 and 0 < x < 2, we get the plot in Figure 
2.9. 
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From Figure 2.9 we can see the extreme value of ^r- increases as £ goes to zero and x goes 
to one. 



Figure 2.9: Power as a function of damping factor and frequency (0.1 <c< 1 and 0 < * < 2) 
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The extreme value of net electric power generated is further illustrated in Figure 2.10 with 
parameters 0.05 < £ < 1 and 0 <x<2. 



Figure 2.10: Power as a function of damping factor and frequency (0.05 <c< 1 and 0 < x < 2). 


2.4 Model Comparison with Experimental Results 

Roundy and Wright [15] use similar mathematical methods to design and test a piezoelec¬ 
trical vibration based generator. Their design uses an 8.44 gram mass and a two-layer sheet 
of PZT-5A with a steel center shim excited by vibrations of 2.5 m s' 2 at 120 Hz- They 
assume the resonance frequency of their generator, co n , matches the driving frequency, co. 
The damping ratio, £, is measured by applying an impulse to the system, and then measur¬ 
ing the open circuit voltage output. The resulting damped harmonic oscillation is used to 
calculate a damping ratio near 0.015 [15]. Figure 2.11 shows the piezoelectric generator 


21 





































designed by Roundy and Wright. 



Figure 2.11: Piezoelectric generator prototype used by Roundy and Wright [15]. 


With these parameter values, we can compute the power generated using our ODE model. 
All red text below are values taken from the Roundy and Wright generator model. 


(0 — (on — {2 k)\20Hz 


second 


Amplitude 


co 

x— — = 1 
(O n 


^ d d 
dc 2 moon 


= 0.015 


m k, 8.44 grams = 8.44 x 10 i kg 
2.5 m 

acceleration = 


>o = 


acceleration 


ftp 


second 2 

2.5 m 
second 2 

( 240 K y 
\ second ) 


= 3.518 x 10' 


in 


(2.36) 

(2.37) 

(2.38) 

(2.39) 

(2.40) 

(2.41) 
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C - mOo 2 ®n = (8-44 x l(T 3 kg)(0.015)(3.518 x 10 ~ 6 m) 2 ( 24Q M (2.42) 

\second J 


C = 


1.0494 x I0~°kg ■ tn 2 


(2.43) 


second 3 

Substituting values given from Roundy and Wright and values calculated above into Equa¬ 
tion (2.29), we get: 


n i) = 


1.0494 x I0~ 6 kg ■ m 2 
second 3 


1 


1.2 x \0~ 3 kg-m 2 


(2(-015 )Y 


second 3 


(2.44) 


The corresponding result from Roundy and Wright is 

3.75 x 10 ~ 4 kg-m 2 


P( 1) 


second 3 


(2.45) 


Roundy and Wright’s experiment results in power 68.75% less than our mathematical re¬ 
sults. The disparity between our ODE model calculation and Roundy and Wright’s exper¬ 
imental results can be explained through the fact our model is purely mechanically based 
and does not take into account loss of energy through the transfer of mechanical energy 
to electrical energy. Despite the power disparity, Roundy and Wright’s results support our 
ODE model. 


2.5 One-Dimensional Power Generation Model Conclu¬ 
sion 

Our 1-D model shows potential to serve as a foundation to optimize power generated by 

p( x \ 

a piezoelectric generator. In order to optimize power generated, -^4, we should minimize 
the damping factor, £, and choose the optimizing frequency ratio, x = Keep in mind 
that as £ approaches zero the optimal frequency ratio approaches one and the maximum 
power generated drastically increases. 
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CHAPTER 3: 
Single-Mode Model: 
Derivation of the Six Variable Model 


There are many studies on the transduction of electrical energy from a single piezoelectric 
transducer (see [18] and references therein) and the general conclusion is that piezoelectric 
energy harvesting approaches have considerable promise for applications involving small 
amounts of power [19]. To boost the power output of such system, one proposal by Scruggs 
[20] is to use multiple transducer patches and channel energy generated by some of the 
patches into other patches in order to excite additional vibration modes. An example of 
this type of energy harvesting system is depicted in Figure 3.1. The parameters describing 
the size of a piezoelectric patch take the typical values l\ « h ~ I 3 — 100 mm, w = 25 mm, 
t p = 0.25 mm and the thickness of the beam is //, = 0.25 mm [21]. 



Figure 3.1: Energy harvesting bimorph cantilever with distributed piezoelectric transducers (gray) 
bonded to a substrate (white). Adapted from Dutoit, Wardle and Kim [16]. 
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The general equations for linear piezoelectricity are shown in Figure 3.2. 


Gauss' law : 

[A 

Z>|>)= d 2 


cD 


cx. 


= P: 


a vector of electrical displacements (charge area) 


LA. 

p. charge density 

Mechanical energy balance: yC 

u t : displacement. T.. : stress 

Constitutive relation [2]: 


c u. 
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c n 

C 12 

C 13 



A 

t 2 

= 


“ e 31 

C 12 


C 13 



A 

t 3 



-e-3 

c 13 

C 13 

c 33 



A 

A 



~ e l5 




C 55 


A 

A 


~ e lS 






c 55 

A 

1 — 

Ok 1 

1_ 


. o 
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= 
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A 
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.A. 
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material stress (force area). 


E(x) = 


A 

A 

A 


electrical field in the material (volts meter) 


e.. : coupling coefficients £.. : : dielectric constants c. : modulus of elasticity 


Figure 3.2: General equations for linear piezoelectricity [20]. 
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The constitutive relation in Figure 3.2 is for material that is electrically and mechanically 
isotropic in the 1 and 2 directions, while the 3-direction is the poling direction. The first 
three rows of the constitutive relation relate the electrical displacement to the electrical field 
and to the strain in the material, and the remaining six rows relate the mechanical stress to 
the electric field and to the strain. Here the stress, T(x), and strain, S(x), are represented in 
Voigt notation or Voigt form. Referring to the Figure 3.1 above, 1 corresponds to x and 3 
corresponds to y. 

Equations for a patch on a beam are largely based on Euler-Bernoulli (or Timoshenko) 
beam models. Piezoelectric layers or patches are bonded to the both sides of the beam. The 
beam is assumed to be in pure bending; all other deformations are considered negligible. 
Here we consider a “bimorph” configuration. That is, in the undisturbed state the beam lies 
along the x-axis in a <x < b,— % <y < — y < z < f, with equal patches on its upper 

and lower surfaces y = each of thickness t p and covering the entire upper and lower 
surfaces. It is further assumed that the patches are bonded perfectly to the substructure 
and the electrodes cover their entire upper and lower surface areas. The beam is clamped, 
at x = a, on a base that is subject to motion in the y-direction; it is free at x = b. The 
polar direction on the top patch is in the y-direction and on the bottom patch it is reversed 
so that the voltages on the two patches add. Due to this asymmetry, we can compute the 
piezoelectric effect on the upper patch alone and double the result. 

Adopting the standard assumptions of Euler-Bernoulli beam theory, a balance of forces and 
moments can be combined to yield 

Mxx = wi (®tt T htt ) • (3.1) 

where M(x, t ) is the internal moment generated by mechanical and electrical strain, m is the 
mass per unit length of the composite beam, co is the relative vertical deflection along the 
beam from the base, and h is the absolute motion of the base/host structure (thus w + h is 
the absolute displacement of the beam). The Euler-Bernoulli bending deformation model 
relates the stress to the curvature of the beam: 


Si = -y(O xx 


(3.2) 
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The constitutive relations yield T\ = cnSi =F^3t£3, D 3 = ±^31^1 + ^33£3 (the sign of 
the coupling terms may vary from patch to patch depending on whether they are poled 
in the positive y or negative y direction). Additionally, it is assumed that the electrical 
field is constant throughout the thickness of the patches; therefore, £3 = ±p with the sign 
depending on the poling direction, where v p is the terminal voltage on the top of the top 
patch with thickness t p . 

The relation between moment and stress is: 

[b 
2 

M = — w I T\ydy—w 

_'± 

2 

M = c 1 \w ~~ (O xx - e 3 iw — + tp [H(x-a)-H (x - b)} v p ■ 2, (3.4) 

where w is the width of both the beam and the piezo patches and H denotes the Heaviside 
step function. 

When the electric field is normal to the beam axis and uniform in the patch over a < x < 
b and zero outside this interval, the piezoelectric effect contribution to the displacement 
equation enters as derivatives of delta functions at the ends of each patch [ 22 ]: 

m(O tt + (O xxxx + e 3 ] w(t h + t p ) [ 8 '{x-b) - 8 ' (x-a)] v p = -mh tt . (3.5) 



The boundary conditions for equation (3.5) are Co(a,t ) = CQ x (a,t ) = 0. This reflects the fact 
that the beam is clamped at x — a, but at the free end, x = b, we have M = M x — 0. The 
size of the charge can be found by integrating the Gauss’ law over the surface area of the 
electrodes [22]: 


q(t)= // D 3 dA - D 3 dA = *? 3 i w(t p + t s ) [eoj 

JJupper JJlower 


2C33 


w(b — a)v p . (3.6) 


Equations (3.5) and (3.6) provide a complete system of equations for the vibration energy 
harvester. The most common method of solving the system (3.5-3.6) is to assume that the 
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beam deflection can be represented by an eigenfunction expansion 

oo 

ft) (x.t) = Y, fyi ( x ) Hi (0> (3.7) 

i= 1 

where (j)i(x) is the i-th mode shape function and «,•(/) is the i-th modal displacement. The 
mode shape functions must satisfy the boundary conditions and follow the following or¬ 
thogonality relationships: 

b 

I (j>i (x) <j>j (x) dx = 8 ij. (3.8) 

a 

Substituting (3.7) into (3.5), multiplying both sides by (j) k (x), and integrating v from a to b 
yields for each k, we arrive at 


W (0 + [wscm] u k +A k Vp = -m [' Yk} h" (t), (3.9) 


where Equation (3.10) is the modal short-circuit frequency, Equation (3.11) is the modal 
electromechanical coupling coefficient, and Equation (3.12) is the modal influence coeffi¬ 
cient of the base excitation. 


w SC,k 


V 


Cii wt b 3 
12 m(b — a) 


(3.10) 


where X k is the eigenvalue corresponding to ^.(v); for cantilevered configurations, these 
numbers are the solutions to cos(A) cosh(A) = 1 


A k = c 3 iw {t b + tp) [$k {x-b)~ Qk {x - a )], (3.11) 

b 

Yk = J Qk(x)dx. (3.12) 

a 

At this point, it is common to add in a modal damping term 2wsc.k^k u k (0 to the left side 
of Equation (3.9) to account for all proportional damping effects. Four possible damping 
mechanisms, one external and three internal, and various combinations of these mecha¬ 
nisms have been considered by Banks and Inman [23]. They are viscous air damping, strain 
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rate damping or Kelvin-Voigt damping, spatial hysteresis, and time hysteresis, respectively. 
Their study shows that air damping plays a more significant role in lower modes whereas 
internal damping plays a more important role for higher modes. 

Equations (3.9) and (3.6) provide a complete system of equations for the vibration en¬ 
ergy harvester. If the beam is excited with a driving frequency that is close to one of 
its natural frequencies, then the corresponding mode of that frequency will dominate 
the motion of the beam. In this situation, it is reasonable to make the approximation 

oo 

(o{x,t) — £ <j)j(x)ui(t ) ~ (f>i (x)u\ (t). After we drop the ”1” subscript, this simple ap- 

(= i 

proximation reduces Equation 3.9 to 

u" (t) + 2w S c^u (; t) + [wsc 2 ] u +Av p = -m [y] h" (t). (3.13) 


Here m is the mass/unit length, y is the modal influence coefficient of the base excitation, 
£ is the modal damping ratio, wsc is the modal short-circuit natural frequency, and A is the 
electromechanical coupling coefficient (effective piezoelectric constant). Equation (3.6) is 
simplified to 


q(t) =e 31 w(t p + t s ) [<M 


b 

a 



a)v p . 


(3.14) 


Differentiating Equation (3.14) yields 


— Au -|-— w (b — a) Vp = —i, 


(3.15) 


where i represents the current. We can focus on the single-mode model (Equations (3.13) 
and (3.15)) for now, which assumes that the beam is vibrating near its fundamental natural 
frequency and, consequently, the motion of the beam is described by one modal coordinate 
(typically, the fundamental mode). We need to solve the system and find the maximum 
power points of energy harvesting over a range of base excitation frequencies, patch thick¬ 
ness and length, load resistances, etc. The results will allow a designer to choose the 
optimal resonant frequency and patch dimensions to maximize the power harvested. 


In the single-mode model, i.e., Equations (3.13) and (3.15), the two unknowns are the 
modal displacement u(t) and the terminal voltage v p , and the two inputs are the base motion 
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h(t) and the current i. The beam tip deflection, which is a predominant measure, can 
be computed by the simple relation c d(L,t ) = (p{L)u{t ), where <j){x) is the fundamental 
mode shape function. The coefficients of Equations (3.13) and (3.15) can be calculated 
using knowledge of the materials and the layout of the beam, or they may be estimated 
experimentally. 


The most commonly analyzed scenario for vibration energy harvesting is steady state base 
excitation. At steady state, we suppose that the beam and base motions are of the following 
forms, respectively 

u(t) = Usm{(O base t), (3.16) 

h it) H sin i(Ob ase t (ftbase ) > (3-17) 

where U is the amplitude of the modal displacement, H is the amplitude of the base motion, 
(Obase is the base motion driving frequency, and (p base is the phase lead of the base motion 
relative to the beam motion. If we further consider a simple resistive load where i = ^ and 
R is the load resistance, then Equations (3.13) and (3.15) are simplified to a linear system. 
Assuming the voltage signal takes the form v p = V p sin((%,„./ — d vo i tage ) and substituting 
these expressions into Equation (3.15), we have 


AU (Obase COS ifObaset) 3“ CqV pG)base COS ( (O b ase^ dvoltage ) 


Vp sin {(Obaset dvoltage ) 

R 


(3.18) 

where Co = f 33 w {b — a). With the help of trigonometric identities cos (x — y) = cosvcosy + 

T p 

sin^siny and sin(x —y) = sin.rcosy — cos.rsiny we can rewrite Equation (3.18) as 


AU (Dbase^OS ( (^baset ) T CoVpObase [cOS (( Qbaset ) COS d voltage T SIA i^baset) sin0 vo /f fl g e J 
= [sin (( Qbaset ) COS dvoitage COS ((Qbaset ) sin d V oltage\ • 

(3.19) 

Matching the coefficients of cos(ftWseO and sin(ft>/ WSY ./) respectively, we obtain 

V p 

AU (Obase + CqV p (Obase COS O vo ltage = sin dvoitage- (3.20) 

Solving for cos d vo i tage , we obtain 


COS dvoitage — CoCObaseR SU1 d V oltage- 


(3.21) 
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Substituting Equation (3.21) into Equation (3.20), we get 


AU(libase "3 Q)VpO)l,ase ( C()Q}f )ase R Qyoltage) ^ sin0 voltage■ (3.22) 


Solving for sin 9 vo i tage , we obtain 


'voltage 


AU (ObaseR 

V p (l+C 0 2 (O base 2 R 2 )' 


Substituting Equation (3.23) into Equation (3.21) yields 


(3.23) 


svoltage — 


CpAU (Phase 2 R 2 
V p (1+Cq 2 (Obase 2 R 2 ) 


Now substituting Equations (3.23) and (3.24) into Equation (3.13) 

U (Obase sin ( (Obaset') 3“ 2w SC C U(Obase COS ( (ObaseR 
3“ [wsC ] U sin {(Obaset) 3“ A Vp si n [(Obaset @voltage ) 
= 171 [y] H(Obase sin ( (Obaset tybase) i 


(3.24) 


(3.25) 


which leads to 

U(Obase sin ( (Obaset ) 3“ 2vV5(7 C,U(Obase^OS (( Obaset ) 3” [wsC - ] U sin (( Obaset ) 
3 ~AVp [sin ( (Obaset COS Oyoltage COS ((Obaset') Sin ^voltage] 

/U [y] H(Obase" [sin (( Obaset) COS (pbase COS (( Obaset ) sin (phase] ■ 

Matching the coefficients of cos ((Obaset) and sin ((Obaset), respectively, we get 


(3.26) 


U(Obase 3“ ['ESC ] U -\-AVp COS ^voltage m [T] H(Obase COS (pbasei 


(3.27) 


2wscCU(Obase AV p sill 0 vo hage — /7i [y] H(Obase sit Mpba 

Substituting Equation (3.21) into Equation (3.27), we obtain 


(3.28) 


U (Oba 


CqAU (Obase" R" 


: . r 21 rr 11/ ™ u, ba se “ 

+ [wsc u + av p^T (1 + c„ 2 ffl,„„ 2 R 2 ) 


= m [y] H (Obase 2 COS (p ha 


(3.29) 
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Solving for cos (j)b ase , we find 


2 | [,,, 21 | C§A (Ojy ase R 

-«W+K| + ( i5S T ? j 

COS (phase = - r^T} - 9 - U. 

111 \y\ HfOhase^ 

Substituting Equation (3.23) into Equation (3.28), we obtain 


(3.30) 


AU (OhaseR 

fp V p (1 + Co 2 (O base 2 R 2 ) 


2 w sc C U (o base +AV P w ^ = ~ m [7] H(O base 2 sin (p ha 


(3.31) 


Solving for sin^ aie , we find 


2 wsc^ 0 ) h 

ase + 


(MbaseR 


sin (p ba 


(l+C() 2 (°base 2 R 2 ) 


m [ 7 ] H (O ba 


U. 


(3.32) 


Using Equations (3.30) and (3.32) and the trigonometric identity cos 2 ((p base ) + sin 2 ((p base ) = 
1, we find 


m, 2 4- L, 21 4- Q)A 2 m base 2 R 2 
+ \_WSC \ + 

/77 [ 7 ] H 0) base 2 


1 2 


u 2 + 


2w5cC®/ 7 + 


^ ®baseR 


( l+C'o^CUiflse 2 ^ 2 ) 


777 [ 7 ] H 


f/ 2 = 1. 


(3.33) 


Solving for U leads to 


jj ///[y]//t0/, fJV( . (l ■ (\) iO] nlSf , 7?“) 

\/( ( ’ [^5C _ ] ) ( 1 SC()~ (Q hase ~R- ) l C()A - (Ob ase "R- ) +(2 W$C^> (Obase ( t ®base~ ~ 

(3.34) 

Using Equations (3.23) and (3.24) and the trigonometric identity cos 2 (6 vo i tage ) + sin 2 (0 vo i tage ) 
1, we obtain 


AU 0J b aseR 

2 

1 

CqAU (Ohase 2 R 2 

Vp (l+Co 2 (d base 2 R 2 ) 

\ 

V p (l +C{ 2 0i base 2 R 2 ) 


Solving for V p , we find 


Vp — AU (O base R. 


(3.35) 


(3.36) 
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Substituting Equation (3.34) into Equation (3.36), we obtain 




ARiii ylH (O b ase (l+Co C0fr ase ~R 


(( ®base~~\~ ^5C“] ) ( 1 +Q) ®base CqA ®base 2 R ) few sc C ®base ( 1 “t“Q) ® base ~R ) ~t~A 2 ®base^) 


(3.37) 


The average power dissipated by the load resistor is given by 


V 2 

p = -A 

2 R 


(3.38) 


Substituting Equation (3.37) into Equation (3.38), we obtain an expression for power: 

P_ _ A 2 Rm 2 [y] 2 H 2 a w 6 (1 +Co 2 co ba J R 2 ) 2 _ 

“ [( ( ^base 2 I [ lv \S'C^] ) ( 1+Q) (Qbase^'R-') A~CoA 2 G)i ) ase~R‘‘') ®base( 1 CQ base ~ R~^) -\-A~ 

(3.39) 

Equation (3.39) provides an explicit and complicated form for power generation, which 
involves nine parameters. We will study the effects of these parameters in the following 
chapters. 
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CHAPTER 4: 

Optimization Methodology 


4.1 Optimization Methodology Overview 

We optimize piezoelectric power generation by analyzing two power generation models 
using two sampling approaches and one iterative approach. The goal of our methodology 
is to discover the optimal values of each variable in order to maximize power generation. 
In the following chapter we compare the results of each optimization method and analyze 
the robustness of our results. 


4.2 Power Generation Models 

We optimize piezoelectric power generation by optimizing the mathematical models de¬ 
scribed in Chapters 2 and 3. We attempt to find the optimal combination of variable values 
in order to maximize power. 


The first model is an adaptation of Equation (2.28). It consists of three variables, ft), ft),,, 
and £. We will not consider m or Yq as they occur only in the numerator of Equation (2.28) 
and obviously maximize power by taking on the largest value possible: 


P C(f ) 3 *> 3 
mY ° 2 (l-g) 2 + (2Q) 2 


(4.1) 


The second model is an adaptation of Equation (3.39). It consists of six variables. We 
will not consider m, y, or H as they occur only in the numerator of Equation (3.39) and 
obviously maximize power by taking on the largest value possible: 

p __ A 2 Rco base 6 { 1 +Cq 2 co base 2 R 2 ) 2 _ 

(tnyH) 2 [(( (O base ~-\-Wsc~)(l^~Co~ (0 ba se~R~)~^C()A-CQ base ~R-)--\-(2wsctD0) base (l-\-CQ~ CQ baS e~R~)-\-A-(Q b aseR)2] 

(4.2) 

Equations (4.1) and (4.2) lie in the core of our analysis. 
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4.3 Variable Ranges 

In order to optimize maximum power generated we must constrain our variables. The goal 
of the optimization is to find the optimum combination of variables in order to maximize 
power and lead to the development of more efficient piezoelectric generators. Keeping in 
mind the practical application of our optimization, we select the ranges of each variable in 
such a way that engineers will be able to select materials and conditions within the limits of 
each variable to create a physical model based on our results. Figures 4.1 and 4.2 display 
the ranges for variables in Equations (4.1) and (4.2), respectively. 

Table 4.1: Variable ranges for the three-variable power generation model of Equation (4.1). 


Variable 

Variable Symbol 

Range 

Frequency of Excitation 

ft) 

[\20ks~\360ks~ 1 ] 

Natural Frequency 

(O n 

[120^" 1 ,360^- 1 ] 

Modal Damping Ratio 

c 

[0.005,0.02] 


In Roundy and Wright’s experiment [15], they use the values co = (O n = 240 ns~ l and £ = 
0.015. We choose the range of ft), ft),, and £ by varying Roundy and Wright’s experimental 
values. 


Table 4.2: Variable ranges for the six-variable power generation model of Equation (4.2). 


Variable 

Variable Symbol 

Range 

Electromechanical Coupling Coefficient 

A 

[0.01,0.99] 

Load Resistance 

R 

[0.02D, 100D] 

Base Motion Driving Frequency 

(&base 

[1207T^” 1 ,3607T^- 1 ] 

Net Clamped Cap. of Piezoelectric Material 

Co 

[1 x 10 8 F, 1 x 10- 6 F] 

Modal Short Circuit Net Frequency 

Wsc 

[120^” 1 ,360^- 1 ] 

Modal Damping Ratio 

c 

[0.005,0.02] 


We choose the range of (Of, ase and £ by varying the experimental values found in Roundy 
and Wright’s paper [15]. Likewise, we choose the range of Co by varying the value cited 
in Fleming and Behren’s paper [24]. We choose to use the same range for wsc as the 
other frequency, OUbase- For the coefficient, A, we use [0.01,0.99]. Finally, the range of 
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R varies from the resistance of a copper wire to about half the resistance of an operating 
incandescent light bulb. 


4.4 m k Factorial Sampling 

One method to find an optimal solution for power generation using our models is to evaluate 
the combination of all independent variables, which we will interchangeably refer to as 
factors in this chapter, at different levels. Levels are the variety of values a factor may 
assume in a simulation. This methodology can be described as m k factorial design, where 
m is the number of levels per factor and k is the number of factors. In an nr factorial 
design the number of design points is equal to m k . Design points are the combination of 
factor levels to be evaluated in the simulation. The larger the value of m for an nr factorial 
design, the better its space-filling properties [25]. However, an increase of levels may 
drastically increase the number of design points required to evaluate in a simulation. 


4.5 Nearly Orthogonal and Balanced Latin Hypercube 
Sampling 

An alternate and more efficient sampling methodology is to use space-filling designs to 
thoughtfully choose the specific levels of each factor to evaluate. Latin hypercube designs 
provide a flexible way of constructing efficient designs for qualitative factors. They have 
some of the space-filling properties of m k factorial design, but require orders of magnitude 
less sampling [25]. 

K.Q. Ye [1998] describes Latin hypercubes as follows: 

An nx d Latin hypercube consists of d permutations of S n = {1,2, ...n}. 

A Latin hypercube design takes row vectors as the experimental points in 
a J-dimensional space. One-dimensional projections of a Latin hypercube 
design are evenly spaced and have the lowest possible discrepancy. [26] 

The significant benefits of the Latin hypercube sampling used in this paper is the samplings’ 
orthogonality and space-filling property. 
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4.6 Sampling Metric: Orthogonality 

The correlation between two vectors u — [u\,u.2, •••««] and v = [vi,V2,...v„] is defined the 
value of the quotient 

L(Vj-v)(Ui~u) 

a/L(m/-m) 2 L(v ; --v) 2 ’ 

where u = ^ and v — If w and v are uncorrelated they are said to be orthogonal 
and their correlation value is zero. Alternatively, two completely correlated vectors have 
a correlation value of one. Orthogonal Latin hypercubes create column vectors with zero 
correlation [26]. 

Orthogonality in variable sampling is particularly valuable in regression analysis and sim¬ 
ulations where one desires independent variables in a regression model to be orthogonal to 
each other in order to minimize variable correlation [27], [28]. Additionally, the orthog¬ 
onal Latin hypercube design ensures independence of estimates of linear effects of each 
variable [26]. 

No matter how many levels or factors, nr factorial samplings result in completely orthog¬ 
onal design points. In other words, the design columns in factorial design are uncorrelated. 

We use a select set of columns from a nearly orthogonal and balanced (NOB) design, 
developed by H. Vieira Jr. [29], in order to reduce the number of design points necessary to 
efficiently find the optimal solution to our power generation problem. After we input our 
model variables, with associated variable ranges, Vieira’s NOB design spreadsheet outputs 
512 nearly orthogonal and balanced design points [30]. According to Vieira, 

Nearly orthogonal means that the maximum absolute pairwise correlation 
between and two design columns is minimal. Nearly balanced means that 
for any single factor column, the number of occurrences of each distinct 
factor level is nearly equal. 

Vieira and Sanchez call a design nearly orthogonal if the absolute pairwise correlation 
between any two factors is less than 0.05. Likewise, they call a design nearly balanced if 
the number of objects in each level of each factor differs from the ideal by no more than 
20% [25], 
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The balance is a concern when discrete-valued factors with limited numbers of levels are 
present because rounding the levels appropriate for continuous-valued factors may induce 
correlations that are larger than desired. Because the designs we use in this paper involve 
only continuous-valued factors, hence the subset of columns is a Latin hypercube, we will 
refer to these as nearly orthogonal and balanced Latin hypercube (NOBLH) designs. 

Below are the correlation between factors for the nearly orthogonal and balanced Latin 
hypercube sampling, consisting of the 512 design points generated by Vieira’s NOB Design 
Spreadsheet. Refer to Appendix, Section A.l and Section A.2 to see the MATLAB code 
we used to generate orthogonality values cited in Tables 4.3 and 4.4. 


Table 4.3: Orthogonality of three-variable NOB Latin hypercube design. 



ft) 


C 

ft) 


-0.0011 

-0.005 


-0.0011 


-0.0038 

c 

-0.005 

-0.0038 



Table 4.4: Orthogonality of six-variable NOB Latin hypercube design. 



A 

R 

®base 

Co 

Wsc 

c 

A 


-0.00106 

-0.0036 

-0.00293 

0.001528 

0.001367 

R 

-0.00106 


-0.0038 

-0.00267 

-0.00195 

-0.00169 

®base 

-0.00359 

-0.00378 


-0.0017 

0.002429 

0.001984 

Co 

-0.00293 

-0.00267 

-0.0017 


-0.00253 

-0.00179 

Wsc 

0.001528 

-0.00195 

0.00243 

-0.00253 


0.006878 

c 

0.001367 

-0.00169 

0.00198 

-0.00179 

0.006878 



Recall that perfect orthogonality is defined as having a correlation value of zero. Notice the 
minimum absolute pairwise correlation between each design column. Using our variable 
inputs, Vieira’s NOBLH design produces close to orthogonal design columns for both our 
three- and six-variable models. 


39 



4.7 Sampling Metric: Space-Filling Property 


It is desirable to select sample design points with good space-filling properties, in other 
words, design points distributed throughout the entire experimental region. For this paper 
we rate the space-filling of the samplings using the Euclidean Maximin (Mm) distance. 

For a given design, define a distance list d = (d\,d. 2 , ...,d[ n ( n _t)]/2)> where 
the elements of d are the Euclidean inter-site distance of the n design points, 
ordered from smallest to largest. The Euclidean Mm distance is defined as 
d\, where a larger value is better. A large value of d\ means that no two 
points are close to (within d\ of) each other. [27] 

We use the Equation (4.4) to compute the Euclidean Mm distance: 

k 

d( s,t) = , | Sj-fj\ 2 (4.4) 

V =1 

where the design space is an n x k matrix and s and t are any two design points [31]. 

When comparing the Euclidean Mm distance of two samples with equal number of design 
points, it is possible to depend solely on Equation (4.4). However, mr factorial and NOBLH 
designs result in differing numbers of design points. NOBLH designs result in 512 design 
points, while the number of m k factorial design points is a function of levels (m) and factors 
(k). For this Euclidean Mm distance calculation we use a three-variable factorial design 
with 10 levels and a six-variable factorial design with five levels, resulting in 1000 and 
15,625 design points, respectively. 

In order to compare Euclidean Mm distance of NOBLH and mr factorial designs, we first 
normalize the design columns of each design with values from the m k factorial design. 
After column normalization, it is possible to use Equation (4.4) to compute the Euclidean 
Mm distance and compare the results. Refer to Appendix, Section A.3 for our Euclidean 
Maximin Distance Code. 

The three-variable and six-variable NOBLH designs result in Euclidean Mm distances of 
0.018 and 0.1947, respectively. The three- and six-variable m k factorial designs result 


40 



in Euclidean Mm distances of 0.1111 and 0.25, respectively. In both the three- and six- 
variable models, the m k factorial designs provide space-filling properties optimal to the 
NOBLH designs. 

Of course, another way of judging the space-filling property of a design is by observing 
the location of the design points in reference to the factor ranges. Unfortunately, this vi¬ 
sual technique is only possible in three dimensions or below. See Figure 4.1 for a visual 
depiction of the 1000 design points evaluated using 10 levels (m) for the three-factor (k) 
model. 



400 iioO 


Figure 4.1: m k factorial sampling. 
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Figure 4.2 visually depicts the 512 nearly orthogonal and balanced Latin hypercube design 
points evaluated for the three-variable model. 



Figure 4.2: Nearly orthogonal and balanced Latin hypercube sampling. 


As the reader will notice, the NOBLH sampling method samples extremely well from the 
interior of the design space and does a poor job sampling from the exterior of the design 
space. The m k factorial sampling method, on the other hand, samples uniformly across the 
design space. 
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4.8 Sampling Metric: Efficiency 

After evaluating all design points in our model, we determine the maximum power gener¬ 
ated and determine an optimal combination of factor levels. This nr factorial method to 
evaluate the six-variable model is easily accomplished using the MATLAB code referenced 
in Appendix, Section A.6. 

The six-variable power generation model requires 49 floating point operations (flops) in 
order to compute one design point. The MATLAB code in Appendix, Section A.6 evaluates 
the power generated for all combinations of the six factors at a varying number of levels. 
This equates to m 6 design points. Therefore the efficiency of the algorithm with the given 
number of factors and levels is 49 x m 6 flops. This model only consists of six factors and 
therefore the number of design points required does not change as dramatically with an 
increase of levels as a model with more factors. For this paper, we use as many as 15 levels 
per factor. When evaluating the model with 15 levels per factor, 15 6 = 11,390,625 design 
points, and 49 x 11,390,625 ~ 558 million flops, our Core i7 equipped computer can carry 
out the algorithm in just over two minutes. 

If our model consisted of ten factors, an nr factorial sampling of 15 levels would consist 
of 15 10 ^ 5.77 x 10 11 design points and too many flops for a standard desk top computer 
to execute in a week. This algorithm and m k factorial design are not a wise use of time for 
simulations with a large number of factors, especially if we desire to change the any aspect 
of the factors or levels and reproduce the simulation. 

The real benefit of nearly orthogonal and balanced Latin hypercube design is its efficiency. 
Using Vieira’s NOB design it is possible to evaluate a model of up to 20 k-level discrete 
factors (where k = 2,3... 11 levels) and 100 continuous factors using only 512 design points 
[30]. Many models or simulations contain many more factors and computations. NOB 
designs choose very good samplings to test given the limited number of simulations a 
system can compute in a given time. 

4.9 NOBLH Iterative Method 

An additional optimization method involves multiple iterations of sampling the design 
space and using results to restrict the design space to choose design points. In this pa- 
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per we refer to this method as the NOBLH Iterative method. 

First, we use NOBLH design to choose design points from the original design space set by 
our variable ranges. Next, we analyze the design point that generate the maximum power 
output. By studying the values of each variable in the optimal design point, we are able 
to make educated guesses about how to restrict the variable ranges in order to focus our 
design space. 

With each newly discovered optimal design point we are able to restrict the variable ranges 
and design space until we reach a maximum power or identify an issue that causes a cease 
to iterations. The aspiration of this method is to arrive at similar results to that of the more 
costly m k factorial sampling method while evaluating fewer design points. 

It should be noted, the iterative method we use is not as systematic as other iterative meth¬ 
ods. Through the use of our NOBLH iterative method we wish to gain insights into our 
optimization models and better understand alternative methods to optimization. 
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CHAPTER 5: 

Optimization Results and Analysis 


5.1 Three-Variable Model NOBLH and m k Factorial 
Results 

After sampling our design space using both NOBLH and nr factorial sampling methods we 
use the resulting design points to find the optimal values of the three variables (ft), ft),,, and 
£) in order to maximize power generation. For the nr factorial samplings we use levels 
of 10 and 100. See the Appendix for the Three-Variable Model Optimization Codes for 
factorial sampling (Section A.4) and NOBLH sampling (Section A.5). The results of each 
sampling method are summarized in Table 5.1. 


Table 5.1: Optimization of power in the three-variable model. 


Sampling Method 

Design Points 

Optimal Power 

ft) 

(On 

c 

NOBLH 

512 

2.49 x 10 10 

339.8n 

339.8 K 

0.0122 

m k (m =10) 

1,000 

7.23 x 10 10 

360k 

360 n 

0.005 

m k ( m = 100) 

1,000,000 

7.23 x 10 10 

360 K 

360 K 

0.005 


The reader will notice both nr factorial sampling methods are able to find an optimal 
power close to three times larger than the NOBLH sampling method. The reader will also 
notice both nr factorial samplings determine the optimal values of each factor to be at the 
extreme of the factors’ ranges. The optimal values of ft) and ft),, are the maximum values 
allowed by the range, 3607T, and the optimal value of £ is the minimum value allowed by 
the range, 0.005. As noted in Chapter 4, the NOBLH sampling fails to adequately sample 
the extremes of the design space, while the m k factorial method samples the entire design 
space uniformly. 

This failure of sampling design points in the extremities of the design space accounts for 
only part of the discrepancy of optimal solutions discovered by the sampling techniques. 
Notice ft) = ft),, for both optimal values for the nr factorial samplings and the NOBLH 
sampling. After further investigation into the model, we find the optimal values for the ft) 
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and co n are always equal. 


We can prove this mathematically by analyzing Equation (4.1). In order to maximize the 
model we obviously want to make the numerator as large as possible while making the 
denominator as small as possible. In order to eliminate the first part of the denominator, 
(1 — ^t) 2 , we set co = co n . This reduces the model to the following equation: 


p c ® 3 0)3 

^ 4 ^ “ 4 £' 


(5.1) 


The three-variable model reduces to a two-variable model in which co — co n . In order to 
maximize power, one simply maximizes co and minimizes £. 


After analyzing the design points, we discover that, out of the 512 design points in the 
NOBLH sampling, only four design points set co = CO n , compared to 100 and 1000 points 
for the m =10 and m =100 nr factorial samplings, respectively. 
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Setting co — (O n plays a crucial role in optimizing this model. The m = 10 m k factorial 
design selects 100 design points with co — ( 0 „. Out of the said 100 design points, 16 de¬ 
sign points produce power greater than the maximum power generated by any of the 512 
NOBLH design points. See Figure 5.1 for the sixteen m k factorial design points that out¬ 
perform the NOBLH sampling. 


U’ 



Figure 5.1: m k factorial design points achieving higher power generation than all NOBLH design 
points. 


While optimizing this particular model, rrr factorial sampling proves the best sampling 
method. After analysis, we determine that to maximize power one can set co = (O n and 
minimize £. 
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5.2 Three-Variable Model NOBLH Iterative Results 

After determining the optimal power generation and associated variable values for the 
three-variable model through the m k factorial results, we attempt to replicate the results 
using the NOBLH iterative method. Our initial hopes are the NOBLH iterative method 
would be able to reach similar results while evaluating less design points than the more 
costly nr factorial sampling method. In this particular case, the three level nr factorial 
sampling finds the optimal power using only 1000 design points. Even if the NOBLH iter¬ 
ative method finds the optimal results in two iterations, the method would use 
512 x2 = 1,024 design points to do so. 

The first iteration of the NOBLH iterative method is identical to the NOBLH sampling 
method. Table 5.2 uses the same variable ranges as in Table 4.1 and Table 5.3 arrives at the 
same results as in Table 5.1. 

Table 5.2: Iteration one NOBLH iterative variable ranges. 



ft) 

(fin 

c 

Lower Limit 

120;t 

\20 k 

0.005 

Upper Limit 

360ti 

360 Ti 

0.02 


Table 5.3: Iteration one NO 

BLH iterative results. 

Optimal Power 

ft) 

(fin 

C 

2.49 x 10 10 

339.87T 

339.87: 

0.01222 


After analyzing the results listed in Table 5.3, we are able to restrain the ranges of ft), (O n 
and £, listed in Table 5.4. Using the new variable ranges we produce 512 new design 
points using NOBLH sampling. After evaluating the new design points, we determine the 
resulting optimal power generation and associated variable values, as shown in Table 5.5. 


Table 5.4: Iteration two NOBLH iterative variable ranges. 



ft) 

(fin 

c 

Lower Limit 

318.3zr 

318.3;r 

0.005 

Upper Limit 

360;r 

360;r 

0.015 
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Table 5.5: Iteration two NOBLH iterative results. 


Optimal Power 

ft) 

(O n 

C 

5.38 x 10 10 

349.27T 

349.97T 

0.0055 


Analyzing the results in Table 5.5, we see optimal power increase again along with both 
values of ft) and ft),,. The optimal value of £ continues to decrease. We use this data to 
further tighten all variable ranges for the third iteration, listed in Table 5.6. 

Table 5.6: Iteration three NOBLH iterative variable ranges. 



ft) 

ft),, 

c 

Lower Limit 

340.67T 

340.6;r 

0.005 

Upper Limit 

360;r 

360;r 

0.006 


Table 5.7: Iteration three NOBLH iterative results. 


Optimal Power 

ft) 

(O n 

C 

6.73 x 10 10 

357.17T 

357 An 

0.005135 


Table 5.7 displays an additional increase in optimal power, ft), and ft),,. Additionally, the 
optimal value of £ continues to decrease. We are able to restrain the ranges for all variables 
for the fourth iteration, as listed in Table 5.8. 

Table 5.8: Iteration four NOBLH iterative variable ranges. 



ft) 

(On 

c 

Lower Limit 

353.3k 

353.3 n 

0.005 

Upper Limit 

360 K 

360 K 

0.0052 


Table 5.9: Iteration four NOBLH iterative results. 


Optimal Power 

ft) 

ft),, 

C 

7.12 x 10 10 

359.4;r 

359.5 K 

0.005029 


Table 5.9 displays optimal power close to the results found using m k factorial sampling. 
After analysis, we further restrain variable ranges for the fifth iteration, listed in Table 
5.10. 
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Table 5.10: Iteration five NOBLH iterative variable ranges. 



ft) 

(O n 

c 

Lower Limit 

359;t 

359k 

0.005 

Upper Limit 

360;r 

360 n 

0.0052 


able 5.11: Iteration five NOBLH iterative results. 


Optimal Power 

ft) 

(On 

c 

7.23 x 10 10 

369k 

360 n 

0.005002 


After five iterations we find optimal power and associated variable values (listed in Table 
5.11) equivalent to m k factorial results. Our conservative methodology to iterative range 
restriction forces us to reach our results after analyzing 5x512 = 2,560 design points. The 
three level m k factorial sampling optimization method determines the same results after 
computing only 1000 design points. However, the NOBLH iterative method shows great 
potential when dealing with a large number of factors. 


5.3 Six-Variable Model NOBLH and m k Factorial 

Results 

After sampling our design space using both NOBLH and m k factorial sampling methods, 
we use the resulting design points to find the optimal values of the six variables (A, R, 
(Obase > Co, wsc and £) in order to maximize power generation. For m k factorial samplings 
we use levels of 3, 5, 7, 10, 11, and 12. See the Appendix for the Six-Variable Model 
Optimization Codes for factorial sampling (Section A.6) and NOBLH sampling (Section 
A.7). The results of each sampling method are given in Table 5.12. 
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Table 5.12: Optimization of power in the six-variable model. 


Method 

Design Points 

Optimal Power 

A 

R 

®base 

Co 

Wsc 

£ 

NOBLH 

512 

7.8509 x 10 9 

0.754 

13.129 

342.61tt 

4.1 x 10~ 8 

344.96 n 

0.006 

m k (m = 3) 

729 

1.8082 x 10 10 

0.5 

50.01 

360 ;t 

1 x 10~ 6 

360 n 

0.005 

m k (m = 5) 

15,625 

1.8082 x 10 10 

0.5 

50.01 

360k 

1 x 10~ 6 

360 k 

0.005 

m k (m = 7) 

117,649 

1.8256 x 10 10 

0.337 

100 

360n 

1 x 10- 6 

360 n 

0.005 

m k (m =10) 

1,000,000 

1.8256 x 10 10 

0.337 

100 

360k 

1 x 10~ 6 

360 n 

0.005 

m k (m = 11) 

1,771,561 

1.8168 x 10 10 

0.402 

70.01 

360 % 

1 x 10~ 6 

360 k 

0.005 

m k (m = 12) 

2,985,984 

1.8203 x 10 10 

0.356 

90.91 

360 n 

1 x 10~ 6 

360 n 

0.005 


Like the three-variable model, all m k factorial samplings discover combinations of factor 
values capable of generating power more than twice the value of the NOBLH sampling’s 
optimal solution. Also like the three-variable model, the optimal solution in the m k factorial 
samplings includes four variables at the extremes of their ranges. The optimal values of 
(Obase, Cq, and wsc are the maximums of the their ranges, while the optimal value for £ is 
the minimum of its range. As we noted in the previous section, NOBLH sampling fails to 
sample the extremes of the design space. The simulations using the NOBLH design points 
does not have the opportunity to explore power generation at the extremes of the design 
space. With only 217 more design points than the NOBLH sampling, the three level m k 
factorial sampling creates design points in the extremes of the variable ranges and finds an 
optimal solution more than twice the value of the NOBLH design solution. 

An interesting result of our simulation is that creating a finer grid (increasing the number 
of levels and resulting design points) does not always result in a more optimal solution. 
For example, both the m = 7 and m =10 m k factorial samplings find an optimal solution 
of 1.8256 x 10 10 . However, when the number of levels increases from ten to eleven the 
optimal solution decreases to 1.82168 x 10 10 . Additionally, the reader will notice the max¬ 
imum power found with m =12 sampling is also less than that discovered with m =10 
sampling. This can be explained by the fact both the m = 7 and m = 10 samplings include 
A = 0.337, while the m =11 and m = 12 samplings do not. 

The true optimal value depends on the precise sampling of A and R. Using the assumption 
the optimal power generation value occurs at the maximum allowable values of (Obase , Q), 
and wsc an d minimum allowable value of £, we are able to further refine the true optimal 
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values of A and R. We accomplish this refinement by fixing (Obase, Cq, wsc, and £ at their 
optimal value while sampling 1000 levels each of A and R. This results in refining the 
optimal values of A as 0.3376 and R as 100. With this more precise value for A we arrive at 
a new optimal power generation only 6.73 x 10 4 greater than our previous optimal solution. 

5.4 Six-Variable Model NOBLH Iterative Results 

After determining the optimal power generation and associated variable values for the six- 
variable model through the mr factorial results, we attempt to replicate the results using the 
NOBLH iterative method. Our aspiration is the NOBLH iterative method will be able to 
reach similar results while evaluating less design points than the more costly m k factorial 
sampling method. From Table 5.12, we see using seven levels the mr factorial sampling 
method results with an optimal power of 1.8256 x 10 10 after evaluating 117,649 design 
points. We deem the NOBLH iterative method a success if we reach a optimal power 
close to the seven level mr factorial sampling method results with far fewer design points 
evaluated. 

The first iteration of the NOBLH iterative method is identical to the NOBLH sampling 
method. Table 5.13 uses the same variable ranges as in Table 4.2, and Table 5.14 arrives at 
the same results as in Table 5.12. 


Table 5.13: Iteration one NOBLH iterative variable ranges. 



A 

R 

®base 

Co 

Wsc 

c 

Lower Limit 

0.01 

0.02 

120 n 

1 x 10~ 8 

120 n 

0.005 

Upper Limit 

0.99 

100 

360k 

1 x 10” 6 

360 n 

0.02 


Table 5.14 

: Iteration 

one NOBI 

LH iterative results. 

Optimal Power 

A 

R 

®base 

Co 

Wsc 

C 

7.85 x 10 9 

0.754 

13.129 

342.6 K 

4.1 x 10~ 8 

3457T 

0.005998 


After analyzing the results listed in Table 5.14, we are able to assume optimal power in¬ 
creases when A, R, (Ob ase and wsc take on values closer to their upper most range value. 
We also assume £ and Co taking values closer to their minimum allowable values results in 
an increase in power. From these assumptions we tighten the ranges of all three variables, 
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listed in Table 5.15. Using the new variable ranges we produce 512 new design points using 
NOBLH design. After evaluating the new design points we determine the resulting optimal 
power generation and associated variable values, given in Table 5.16. 


Table 5.15: Iteration two NOBLH iterative variable ranges. 



A 

R 

®base 

Co 

Wsc 

c 

Lower Limit 

0.25 

1 

3407T 

1 x 10~ 8 

3A0k 

0.005 

Upper Limit 

0.99 

100 

360;r 

1 x 10~ 7 

3607T 

0.01 


Table 5.16: Iteration two NOBLH iterative results. 


Optimal Power 

A 

R 

t &bcise 

Co 

WSC 

C 

1.64 x 10 10 

0.812 

13.98 

358.57T 

1.28 x 10 -8 

358.7tt 

0.00533 


After the second iteration we notice an increase in optimal power, (0b ase , and wsc- The 
optimal A and R increase slightly. The optimal £ and Cq reduce slightly. Using these 
findings we further restrain variable ranges, listed in Table 5.17. Our results are shown in 
Table 5.18. 


Table 5.17: 1 

teration three NOBL 

H iterative variable ranges. 


A 

R 

t &base 

Co 

WSC 

c 

Lower Limit 

0.5 

5 

3587T 

1 x 10“ 8 

3587T 

0.005 

Upper Limit 

0.99 

50 

360;t 

5 x 10~ 8 

360k 

0.006 


Ta 

ble 5.18: 

Iteration three NOBLH iterative 

results. 


Optimal Power 

A 

R 

®base 

Co 

Wsc 

c 

1.778 x 10 10 

0.626 

28.24 

358.4;r 

4.1 x 10" 8 

358.2;r 

0.005006 


From the third iteration we notice an increased optimal power along with a slight decrease 
in the optimal values of (Ob ase , wsc, and £. The optimal values of A, R , and Co increase 
slightly. We continue restricting variable ranges, listed in Table 5.19. 
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Table 5.19: 

Iteration four NOB 

LH iterative 

variable 

ranges. 


A 

R 

t &base 

Co 

wsc 

c 

Lower Limit 

0.5 

20 

3587T 

2 x 10' 8 

358tt 

0.005 

Upper Limit 

0.8 

30 

359; t 

5 x 10” 8 

359k 

0.00501 


Tab 

le 5.20: 1 

Iteration four NOBI 

LH iterative results. 

Optimal Power 

A 

R 

(fibase 

Co 

wsc 

c 

1.792 x 10 10 

0.728 

21.31 

358.9;r 

2.9 x 10~ 8 

358.9;r 

0.005 


From the results of the fourth iteration, Table 5.20, we see another increase in optimal 
power. However, knowing the optimal values discovered from m k factorial sampling, we 
have made incorrect assumptions restricting the ranges of A, R, (db ase , Co, and wsc- If 
we further restrict the variable ranges in subsequent iterations we would fail to reach an 
optimal power similar to that found through mr factorial sampling. 

An alternative way to look at the NOBLH iterative results is purely in terms of the number 
of design points evaluated. After four iterations, we only evaluate 4x512 = 2048 design 
points. The optimal power resulting from the NOBLH iterative method (1.792 x 10 10 ) 
only represents a 1.84% loss in power from the seven level mr factorial sampling method 
(1.8256 x 10 10 ). Less than a two percent loss in optimal power seems minor when you 
consider the seven level m k factorial sampling method requires the evaluation of 117,649 
design points. Again, the NOBLH iterative method shows potential when optimizing mod¬ 
els with a large number of variables. 

5.5 Robustness Analysis 

The sampling methods provides the means to find the optimal combination of variable 
values to maximize power generation in our models. After running simulations with the 
sampled design points we find the optimal solution for each model. On paper, we then 
possess the exact metrics to build a piezoelectric generator prototype in order to provide 
optimal power output. However, the precision of the variable values we insist on paper 
are not always possible to replicate in experiments or practical application, even under 
perfect conditions. All variables depending on materials or conditions have associated 
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variability. For example, a resistor may be advertised as providing a certain number of 
ohms of resistance, but in reality, provides the advertised ohms of resistance plus or minus 
some level of discrepancy allowed by the manufacturer. 

The purpose of robust design, as a process of simulation optimization, is to identify the 
“best” solution as one not overly sensitive to small changes in system inputs. In this section 
we will attempt to identify which factors in the design of a piezoelectric generator need 
more precision and which factors will allow more variability while producing as much 
power as possible [32]. 

5.5.1 Three-Variable Model 

Using mathematical models to simulate a real world problem, it is easy to see the variables 
only as their numerical values. In reality the variables respond to specifications of the 
piezoelectric generator and its environmental setting. 

In this simple three-variable model, engineers will be able to control all of the variables 
in experimentation. They will be able to set the frequency of excitation (co) by adjusting 
the environmental conditions affecting the piezoelectric generator. Additionally, they will 
be able to set the natural frequency (co„) and modal damping ratio (£) of the generator by 
adjusting the design of the prototype. 

In the practical application of a piezoelectric generator, certain conditions may not be as 
easy to control. A piezoelectric generator designed to operate using the vibrations of an 
air conditioning vent will not be able to depend on a precise or consistent frequency of 
excitation. Engineers must design the generator to operate within a range of frequency 
while producing the necessary power. 

Earlier in this chapter, we determined that the optimal solution to the three-variable model 
occurred when co = co„. Given the fact the environment determines the frequency of exci¬ 
tation, we fixed co during the robustness analysis while varying co„ and £. 

Single Variable Manipulation 

We have determined the optimal solution of the three-variable model depends on minimiz¬ 
ing £ while maximizing co and setting co = G0 n . Due to this knowledge we evaluate the 


55 



robustness of the optimal solution to the three-variable model in terms of varying £ and 
varying the difference between co and ft),,. 

First, we set co = ft),, = 3607T and vary £ in order to evaluate the power lost while deviating 
from the optimal value of £. See the results listed in Table 5.21. 


Table 5.21: Single variable manipulation of £. 


c 

% change in £ 

Power 

% Loss 

0.005 

baseline 

7.23 x 10 10 

baseline 

0.00505 

1% 

7.162 x 10 10 

0.99% 

0.00525 

5% 

6.889 x 10 10 

4.76% 

0.0055 

10% 

6.576 x 10 10 

9.09% 

0.006 

20% 

6.028 x 10 10 

16.67% 

0.01 

100% 

3.617 x 10 10 

50% 

0.02 

300% 

1.808 x 10 10 

74.75% 


A 1%, 5%, or 10% change in £ results in power losses of less than 10%. A change as 
drastic as doubling the desired value of £ results in only a 50% loss in power. 

Next, we set £ = 0.005, ft) = 360n and vary ft),, to evaluate the resulting power lost. See 
the results in Table 5.22. 
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Table 5.22: Single variab 

le manipulation < 

of co n . 

(On 

% change in co n 

Power 

% Loss 

3607T 

baseline 

7.23 x 10 10 

baseline 

360.367T 

0.1% 

6.948 x 10 10 

3.94% 

359.647T 

-0.1% 

6.962 x 10 10 

3.75% 

360.97T 

0.25% 

5.775 x 10 10 

20.16% 

359.17T 

-0.25% 

5.798 x 10 10 

19.84% 

361.87T 

0.5% 

3.608 x 10 10 

50.12% 

358.27T 

-0.5% 

3.626 x 10 10 

49.87% 

363.6n 

1% 

I .444 x 10 10 

80.04% 

356.47T 

-1% 

1.449 x 10 10 

79.96% 

318k 

5% 

7.154 x 10 8 

99.01% 

342 n 

-5% 

7.160 x 10 8 

99.01% 


The reader will notice the drastic power loss when the separation between co and co n in¬ 
creases. We find that even a 1% discrepancy between co and co n results in a power loss 
greater than a 300% change in £ from its optimal value. Whether the ( 0 n is greater than or 
less than co does not make a difference. 

Another way to visualize the sensitivity of co n and £ is to plot the resulting power in terms 
of optimal power as each variable deviates from its optimal value. See Figure 5.2 below. 
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Figure 5.2: Single variable manipulation and resulting percent optimal power. 


As seen in Figure 5.2, the reduction of power when £ increases from its optimal value is 
linear, with a relatively small slope. On the other hand, when ft)„ increases, resulting in 
an increase in the distance between ft) and ft)„, power drastically reduces within the first 
percentages of change in ft),,. 

Alternatively, as £ decreases from its optimal value, power increases linearly with the same 
small slope. The reader should note, reducing £ from its optimal value sets £’s value out¬ 
side the values originally studied in this paper. As ft),, decreases from its optimal value, 
resulting in an increase in the distance between ft) and ft)„, power drastically reduces sym¬ 
metrically to the power reduction resulting from an increase in ft),, from its optimal value. 

Upon completion of the single variable manipulation robustness analysis of the three- 
variable model, we determine the importance of minimizing the discrepancy between ft) 
and ft),,. £ should be set as low as possible, but the precision of setting £ should be a 
second priority to precisely matching ft) to ft),,. 
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Multiple Variable Manipulation 

Single variable manipulation provides insight into how the departure from the optimal value 
of each variable affects power output when the other variables are fixed at their optimal val¬ 
ues. While single variable manipulation provides valuable analysis, it is seldom realistic to 
assume fixing variable values in practical application. A more realistic analysis comprises 
of treating each variable value as a function of a distribution. The distribution has an upper 
and lower limit to match the possible limits guaranteed by a manufacture, for example, a 
capacitor providing 10±1 ohms of resistance. We then analyze how the variable distri¬ 
butions affect the mean and median of power output. From this information we have the 
ability to make informed recommendations in allowing more or less variability in variable 
distributions. 

We assume that engineers would design the piezoelectric generator in order to optimize 
power output with a specific frequency of excitation in mind. The generator itself cannot 
manipulate the environment and accompanying frequency of excitation. To simulate this 
assumption we set the frequency of excitation, oo, to the optimal value of 360 n. We are 
more interested in finding how the variation in design specifications of the piezoelectric 
generator affects power output. 
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In order to guarantee that variable values stayed within desired limits, we use triangle 
distributions. First, using SAS’s commercial statistical discovery software, JMP, we create 
triangle distributions for both co n and £ with modes equal to their respective optimal values 
(tt),,=36071;, £=0.005) and the limits at ±10% the optimal values. Using JMP, we create 
10,000 sample points from the distributions and fixed value of co and computed 10,000 
power outputs using the three-variable model. The distributions of a> n , £ and the resulting 
power output are illustrated in Figure 5.3. 
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Figure 5.3: (O n ± 10%, £ ± 10%, and resulting power distributions. 


Letting a> n and £ vary 10% from their optimal values results in a mean power output of 
9.976 x 10 9 . This mean power represents an 86% loss from the optimal value. Furthermore, 
the median power output is only 2.09 x 10 9 , a 97% loss from the optimal power value. This 
means 50% of all the design points produce a loss of 97% or more. 
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Next, we use insights derived from the analysis of the single variable manipulation to drive 
our quest for more power output. We see from single variable manipulation the importance 
of setting co = (O n in order to maximize power. In the following distribution we allow co n 
to vary only 5% from its optimal value, while allowing £ to vary 10%. See Figure 5.4 for 
the (0,i, £, and resulting power distributions. 
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Figure 5.4: oo n ±5%, £±10%, and resulting power distributions. 


Tightening the distribution to allow for only a 5% variation of co n while keeping the previ¬ 
ously used £ distribution results in a mean power output of 1.782 x 10 10 , a 75% loss from 
the optimal power value. The resulting median power output for these distributions equals 
7.51 x 10 9 , a 90% power loss. Even if the results are not overwhelming, this restriction of 
variation for (O n leads to an improvement in power. 
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Due to the results we see from slightly tightening the ft),, distribution, we decide to inves¬ 
tigate further by allowing ft),, to vary only 1% from its optimal value while allowing £ to 
continue to vary 10% from optimal. See Figure 5.5 for the ft),,, £, and resulting power 
distributions. 
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Figure 5.5: co„±l%, £±10%, and resulting power distributions. 


Further tightening the distribution to allow for only a 1% variation of ft),, while keeping 
the previously used £ distribution results in a mean power output of 5.131 x 10 10 , a 29% 
loss from the optimal power value. The resulting median power output for these distribu¬ 
tions equals 5.4 x 10 10 , a 25% power loss. This is a drastic improvement in power output. 
Notice the power distribution’s shape changes as well, resulting in an increased likelihood 
of generating power close to the optimal and shifting the median power output above the 
mean. 
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Next, we try to further improve upon our findings by restricting the variation of £ from 
its optimal value. We tighten the distribution of £ to only 5% of its optimal value while 
keeping co, ; ’s distribution to allow an deviation of 1% of its optimal value. See Figure 5.6 
for the (o n , £, and resulting power distributions. 
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Figure 5.6: (0„± 1%, £ ±5%, and resulting power distributions. 


Tightening the distribution to allow for only a 5% variation of £ while keeping the pre¬ 
viously used (o n distribution results in a mean power output of 5.128 x 10 10 , a 29% loss 
from the optimal power value. The resulting median power output for these distributions 
equals the median of the previous distributions at 5.4 x 10 10 , a 25% power loss. We find 
tightening £’s distribution actually leads to a decrease in mean power output. This is due to 
the fact when (0 = (o n as £ approaches zero power increases infinitely. That is, there exists 
a vertical asymptote when £ = 0. With a tighter distribution, the value of C does not have 
the ability to take on smaller values. 

Using multiple variable manipulation robust analysis we find when attempting to optimize 
power it is best to allow as little variation as possible from natural frequency’s (c o „) optimal 
value. This will minimize the discrepancy between (0 and (O n and result in greater power 
output. The restriction of the variation of the modal damping ration (£) from its optimal 
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value is not as important. In fact, we find if the variation of £ would allow a smaller value 
than the desired setting, the model generates even greater power. 


5.5.2 Six-Variable Model 

In the six-variable model, engineers are able to control all of the variables in experimen¬ 
tation. They are able to set the base motion driving frequency ((Obase) by adjusting the 
environmental conditions affecting the piezoelectric generator. Additionally, they are able 
to set the electromechanical coupling coefficient (A), load resistance (R), net clamped ca¬ 
pacity of piezoelectric material (Co), modal short circuit net frequency (wsc) and modal 
damping ratio (£) by adjusting the design of the prototype. 

In the practical application of a piezoelectric generator certain conditions may not be as 
easy to control. A piezoelectric generator designed to operate using the vibrations of an air 
conditioning vent will not be able to depend on a precise or consistent base motion driving 
frequency. Given the fact the environment determines the base motion driving frequency, 
we fix (Obase during the robustness analysis while varying the other five factors. 


Single Variable Manipulation 

Earlier in this chapter, we determined the optimal solutions of the six-variable model in 
order to optimize power. To remind the reader of these values we include Table 5.23. 


Table 5.23: Optimal values for each variable in the six-variable model. 


Max Power 

A 

R 

(Obase 

Co 

Wsc 

c 

1.8256 x 10 10 

0.3376 

100 

360 K 

1 x 10~ 6 

360 K 

0.005 


First, we fix all variables at their optimal values and vary the electromechanical coupling 
coefficient (A). See the results in Table 5.24. 
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Table 5.24: Single variable manipulation of A. 


A 

% change in A 

Power 

% Loss 

0.3376 

baseline 

1.8256 x 10 10 

baseline 

0.30384 

-10% 

1.8057 x 10 10 

1.09% 

0.32072 

-5% 

1.8209 x 10 10 

0.26% 

0.3342 

-1% 

1.8254 x 10 10 

0.01% 

0.3410 

1% 

1.8254 x 10 10 

0.01% 

0.3545 

5% 

1.8211 x 10 10 

0.25% 

0.3714 

10% 

1.8088 x 10 10 

0.92% 


From the data presented, it appears varying A does not drastically decrease power output. 
Varying A by 10% from its optimal value only results in around a 1% loss in power. 

Next, we fix all variables at their optimal values and vary load resistance ( R ). See the results 
in Table 5.25. 


Table 5.25: Single variable manipulation off?. 


R 

% change in R 

Power 

% Loss 

100 

baseline 

1.8256 x 10 10 

baseline 

90 

-10% 

1.8175 x 10 10 

0.44% 

95 

-5% 

1.8228 x 10 10 

0.15% 

99 

-1% 

1.8252 x 10 10 

0.02% 

101 

1% 

1.8259 x 10 10 

-0.02% 

105 

5% 

1.8262 x 10 10 

-0.04% 

110 

10% 

1.8251 x 10 10 

0.03% 


By analyzing the data above we see that varying R changes power by even less than the 
variation of A. The reader will notice negative loss when R takes on the value of 101 and 
105 ohms. This represents a gain of power at those values of R. The reader will also notice 
increasing R past those values to 110 ohms results in a power loss from the baseline. It is 
interesting to note, that had we initially included 105 ohms in the range of load resistance 
to be studied we would find a different optimal value of R. 
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Next, we fix all variables at their optimal values and vary the net clamped capacity of 
piezoelectric material (Co). See the results in Table 5.26. 


Table 5.26: Single variable manipulation of Cq. 


C 0 

% change in Co 

Power 

% Loss 

1 x 10~ 6 

baseline 

1.8256 x 10 10 

baseline 

9 x 10~ 7 

-10% 

1.8223 x 10 10 

0.18% 

9.5 x 10 -7 

-5% 

1.8239 x 10 10 

0.09% 

9.9 x 10 -7 

-1% 

1.8253 x 10 10 

0.02% 

1.01 x 10~ 6 

1% 

1.8259 x 10 10 

-0.02% 

1.05 x 10“ 6 

5% 

1.8274 x 10 10 

-0.1% 

1.1 x 10^ 6 

10% 

1.8292 x 10 10 

-0.2% 


The variation of Co from its optimal value results in small variation of power output. Also, 
as Co increases so does the power. From the results above, we assume if we had included 
a larger value as the maximum limit of Co to study we would derive a larger value of Co as 
its optimal value. 

Next we fix all variables at their optimal values and vary the modal short circuit net fre¬ 
quency (vf 5 c). See the results in Table 5.27. 


Table 5.27: Single variable manipulation of wsc- 


Wsc 

% change in wsc 

Power 

% Loss 

360tz 

baseline 

1.8256 x 10 10 

baseline 

324 k 

-10% 

2.0229 x 10 8 

98.89% 

342;r 

-5% 

7.5408 x 10 8 

95.87% 

356.4;r 

-1% 

9.7642 x 10 9 

46.51% 

363.6 K 

1% 

8.5501 x 10 9 

53.17% 

31%n 

5% 

6.5336 x 10 8 

96.42% 

396 K 

10% 

1.6194 x 10 8 

99.11% 


Variation in wsc results in large power loss. As small as a 1% deviation in modal short 
circuit frequency results in close to a 50% loss in power. A 10% deviation of wsc results in 
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a 99% power loss. The six-variable model is extremely sensitive to variations of wsc from 
its optimal value. 

As in the three-variable model, the difference between two frequencies in the six-variable 
model greatly impact power generation. In this robustness study we fix (Obase at its optimal 
value and vary wsc■ D ue to the fixed (Obase, the relationship between the two frequencies 
and power generation is not obvious to the reader. However, after analyzing Equation 
(4.1), the reader will notice when (Obase — w sc the denominator decreases by eliminating 
( — (Obase 2 + wsc 2 )(l + Cq 2 (Obase 2 R 2 ) ■ The decrease in the denominator of Equation (4.1) 
leads to greater power generation. 

Finally, we fix all variables at their optimal values and vary the modal damping ratio (£). 
See the results in Table 5.28. 


Table 5.28: Single variable manipulation of £. 


c 

% change in £ 

Power 

% Foss 

0.005 

baseline 

1.8256 x 10 10 

baseline 

0.0045 

-10% 

2.0227 x 10 10 

-10.79% 

0.00475 

-5% 

1.9203 x 10 10 

-5.19% 

0.00495 

-1% 

1.8440 x 10 10 

-1.01% 

0.00505 

1% 

1.8075 x 10 10 

0.99% 

0.00525 

5% 

1.7377 x 10 10 

4.82% 

0.0055 

10% 

1.6560 x 10 10 

9.29% 


The variation of £ from its optimal value results in small variation of power output. Also, 
as £ decreases the power increases. From the results above we could assume that if we 
had included a smaller value as the minimum limit of £ to study we would derive a smaller 
value of £ as its optimal value, as in Figure 5.7. 
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Another way to visualize the sensitivity of each variable is to plot the resulting power in 
terms of optimal power as each variable deviates from its optimal value. 



Figure 5.7: Single variable manipulation and resulting percent optimal power. 


Figure 5.7 further illustrates the need to set wsc to its precise optimal value. The effects of 
all other variables are drowned out by wsc’s dominating negative effect. Individually A, R, 
and Co change power output insignificantly even under ±10% deviations from their optimal 
values. £ demonstrates its ability to linearly effect power both negatively and positively. 
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In order to see the effects of the less sensitive variables, see the zoomed in plot in Figure 
5.8. 



Figure 5.8: Zoomed-in single variable manipulation and resulting percent optimal power. 


From Figure 5.8 one will notice, after wsc, the deviations of £ have the largest impact on 
power output. Also notice a change of -10% in £ increases power slightly more than a 
change of 10% in £ decreases power. 

The zoomed-in plot also illustrates the curvature of A’s power curve. The optimal value of 
A would not have changed even if we selected a wider range of A values to study. 

Upon completion of the single variable manipulation of the six-variable model robustness 
analysis, we find modal short circuit net frequency (wsc) proves by far the most sensitive 
variable. A deviation of ±10% wsc from its optimal value results in around a 99% power 
loss. On the opposite end of the spectrum a deviation of ±10% A, R, or Co from their 
optimal values results in a power loss of around 1% or less. Interesting to note, slightly 
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increasing R or Co or decreasing £ results in a power gain. 

Multiple Variable Manipulation 

As previously discussed in the three-variable model robustness analysis, the single variable 
manipulation of the six-variable model provides insight into how the departure from the 
optimal value of each variable affects power output when the other variables are fixed at 
their optimal values. Like the three-variable multiple variable robustness analysis we use 
variable distributions to attempt to identify which variables require more precision and thus 
less variation from their optimal value in order to maximize power. 

We assume that engineers design the piezoelectric generator in order to optimize power out¬ 
put with a specific frequency of excitation in mind. The generator itself cannot manipulate 
the environment and accompanying frequency of excitation. To simulate this assumption 
we set the base motion driving frequency, COb ase , 10 the optimal value of 3607T. We are 
more interested in finding how the variation in design specifications of the piezoelectric 
generator affect power output. 
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In order to guarantee that variable values stay within desired limits we use triangle distribu¬ 
tions. First, using SAS’s commercial statistical discovery software, JMP, we create triangle 
distributions for the five variables in the six-variable model with modes equal to their re¬ 
spective optimal values (A — 0.3376, R — 100, Co = 1 x 10 6 , wsc = 3607T, and £=0.005) 
and the limits at ± 10% the optimal values. Using JMP, we create 10,000 sample points 
from the distributions and fixed value of (0/y ase , and compute 10,000 power outputs using 
the six-variable model. See the distributions of the variables and the resulting power output 
in Figure 5.9. 
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Figure 5.9: A± 10%, R± 10%, Co± 10%, wsc± 10%, £± 10%, and resulting power distributions. 


Letting the five variables vary 10% from their optimal values results in a mean power 
output of 4.4875 x 10 9 . This mean power represents a 75% loss from the optimal value. 


71 










































































































































Furthermore, the median power output is only 1.92 x 10 9 , an 89% loss from the optimal 
power value. This means 50% of all the design points produce a loss of 89% or more of 
the optimal value. 


Next, we use insights derived from the analysis of the single variable manipulation to seek 
more power output. From single variable manipulation we determine wsc to be the most 
sensitive variable in the six-variable power generation model. In the following distribution 
we allow wsc to vary only 5% from its optimal value, while allowing the other four vari¬ 
ables a 10% variation. See the distributions of the variables and the resulting power output 
in Figure 5.10. 
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Figure 5.10: A± 10%, R± 10%, Cq± 10%, wsc±5%, £± 10%, and resulting power distributions. 


Allowing wsc to vary only 5% from its optimal value while allowing the other four vari- 
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ables a 10% variation results in a mean power of 7.6514 x 10 9 , representing a 58% power 
loss. This is a 17% improvement in mean power loss from allowing wsc 10% variability. 
Additionally, the median power output is 5.8 x 10 9 , a 68% power loss. The median power 
loss is a 21% improvement over allowing wsc 10% variability. 


Figure 5.10 illustrates the benefit of constraining the variation of the modal short circuit 
net frequency. Next, we further tighten the variation of wsc to 1% while maintaining the 
other four variables at 10% variation. See the distributions of the variables and the resulting 
power output in Figure 5.11. 
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Figure 5.11: A± 10%, R± 10%, Cq± 10%, wsc± 1%, £± 10%, and resulting power distributions. 


Further tightening the variation of the modal short circuit net frequency results in an even 
greater mean power output of 1.599 x 10 10 , representing only a 12% power loss from the 
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optimal power. Furthermore, the reader will notice the power distribution changing drasti¬ 
cally as well, resulting in a median power greater than the mean power. The median power 
output is 1.7 x 10 10 , a power loss of only 7% from the optimal power. 


We learn the importance of minimizing the variation of the modal short circuit frequency 
from its optimal value. If we cannot guarantee a 1% variation constraint on wsc , is it 
possible to make up power output by tightening the variation of the other four variables? 
We attempt to answer the above question by allowing all five variables a 5% variation to 
compare power output to previous results. See the distributions of the variables and the 
resulting power output in Figure 5.12. 
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Figure 5.12: A ±5%, R± 5%, Co ±5%, wsc±5%, £ ±5%, and resulting power distributions. 


From Figure 5.12 we see a power distribution similar to that of Figure 5.10. In fact, restrict- 
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ing all variables to a 5% variation results in a mean power of 7.656 x 10 9 , representing a 
58% power loss from the optimal power. This mean power is extremely similar to the re¬ 
sults displayed in Figure 5.10. Additionally similar to the results from Figure 5.10, the 
median power output is 5.75 x 10 9 , representing a 69% power loss. We find tightening the 
variation of the variables other than the modal short circuit net frequency results in negligi¬ 
ble power improvement over when we restrict wsc to 5% variation and allow 10% variation 
for the other four variables. 


Our next question asks whether it is possible to generate a best case scenario by restricting 
all variables to 1% variation. Do we see significant improvement over when we restrict wsc 
to 1% variation while allowing the four other variables 10% variation? See the distributions 
of the variables and the resulting power output in Figure 5.13. 


A 





R 


G 



£ 




Power 



Figure 5.13: A ± 1%, R± 1%, Cq =t 1%, wsc ± 1%, C 3= 1%, and resulting power distributions. 
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Comparing the power distribution in Figure 5.13 with the power distribution in Figure 5.11, 
we notice the l ik elihood of producing power close to the optimal increases. The power 
distribution in Figure 5.11, however, displays the possibility of producing power greater 
than the optimal. As we discussed in single variable manipulation, this occurs when £ is 
smaller than its optimal value or Co or R is greater than its optimal value. 

Restraining all variables to a 1% variation results in a mean power of 1.601 x 10 10 , repre¬ 
senting a 12% power loss, and a median power of 1.7 x 10 10 , representing a 7% power loss. 
The restriction of all variables to a 1% variation results in a negligible benefit in terms of 
mean and median power compared to restraining modal short circuit net frequency to 1% 
variation while allowing the other four variables 10% variation. 

We see when we tighten the variation of wsc to 1% and allow 10% variation of the other 
variables there is a possibility of generating power greater than the optimal. During single 
variable manipulation we determined decreasing £ from its optimal value has the greatest 
influence in increasing power. We also note decreasing £ 10% increases power slightly 
more than increasing £ 10% decreases power. Are we able to produce a higher mean 
and median power if we allow £ to vary 10% while restricting all other variables to 1% 
variation? See the distributions of the variables and the resulting power output in Figure 
5.14. 
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Figure 5.14: A ± 1%, R ± 1%, Co ± 1%, wsc ± 1%, C T 10%, and resulting power distributions. 


The power distribution in Figure 5.14 looks extremely similar to the power distribution in 
Figure 5.11 when we restrict wsc to 1% variation while allowing 10% variation for all other 
variables. For both power distributions there are possibilities of producing power greater 
than the optimal and a high likelihood of producing power approximate to the optimal. 
Additionally, the power distributions have extremely similar means and medians. The mean 
and median power produced for the distribution above are 1.604 x 10 10 and 1.7 x 10 10 , 
respectively. This accounts for a mean power loss of 12% and a median power loss of 7%. 

Although not an overwhelming front runner, we find that the best mean power results from 
restricting all variables to 1% variation while allowing £ 10% variation. This produces a 
difference of only 3 x 10 7 greater power than restricting all variables to 1% variation. 
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The biggest takeaway from the six-variable robustness analysis is the importance in mini¬ 
mizing the variation of the modal short circuit net frequency (wsc) from its optimal value. 
The restriction of variation for all other variables should be a distant second priority. In 
fact allowing the modal damping ratio (£) to vary from its optimal improves mean power 
production slightly. 
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CHAPTER 6: 
Conclusion 


6.1 Summary of Results 

The military is looking for alternative energy sources to provide deployed bases with elec¬ 
tricity. Diesel fueled generators and their large logistics tails have proven expensive and a 
liability to force protection. The alternative energy focus for now is solar, but solar arrays 
require infrastructure and space. Piezoelectric generators, combined with other renewable 
energy means, could provide a more efficient means of electricity production to fossil fuel 
powered generators. Additionally, new methods of energy storage could eventually facili¬ 
tate alternative energy use in military applications [33]. 

This paper attempted to find the optimal designs involving material parameters in the piezo¬ 
electric generator as well as the generator’s environment in order to maximize electric out¬ 
put. Our hope is that the results of this paper will serve as a step towards the practical 
application of efficient piezoelectric generators and contribute to preserving the combat ca¬ 
pability of the US Military. Through the study of two mathematical models we found the 
optimal values for variables given constrained ranges. We presented the optimal values for 
the three-variable model in Table 5.1 and the optimal values for the six-variable model in 
Table 5.12. 

6.1.1 Three-Variable Model Analysis 

In the detailed analysis of the three-variable model, we study the power as both a function of 
normalized frequency and a function of modal damping ratio and frequency. When we ana¬ 
lyze power as a function of normalized frequency, we find that generating maximum power 
the optimal normalized frequency (x — 77 -) depends on the value of the modal damping 
ratio (£). As £ approaches zero the optimal normalized frequency approaches one. When 
we analyze power as a function of modal damping ratio and normalized frequency, we are 
able to visualize the potential power generation capability by minimizing £ and setting the 
normalized frequency to one. Finally, we use Roundy and Wright’s results [15] to support 
our simple three-variable model. 
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6.1.2 Six-Variable Model Analysis 

Through the analysis of the six-variable model optimization results we find that power in¬ 
creases as the base motion driving frequency (Ct)/, aye ), net clamped capacity of piezoelectric 
material (Co), and modal short circuit net frequency (wsc) increases and modal damping 
ratio (£) decreases. The optimal values for the electromechanical coupling coefficient (A) 
and load resistance ( R ) vary depending on sampling method. There appears to be a corre¬ 
lation in power output and the difference of (Qbase an d wsc- 

6.1.3 Sampling Methods 

In our search for optimal values to maximize power production we use two sampling meth¬ 
ods, nr factorial and nearly orthogonal and balance Latin hypercubes (NOBLH). Although 
the NOBLH design provides nearly orthogonal samplings with fair space-filling properties, 
its failure to effectively sample the extremes of the design space leads to less than optimal 
results. In both of our models, the optimal variable values exists at the extremes of the 
variable ranges in all but one variable (A). m K variable sampling chooses design points 
incorporating the optimal values, while NOBLH does not. 

Had the models required more computational rigor (more factors or computations) the 
NOBLH samplings’ efficiency might have provided insights into optimization where m k 
factorial design would have been too costly. However, with both our three- and six-variable 
models we do not have difficulty evaluating m k factorial samplings even with double digit 
levels (in). 

6.1.4 Iterative Method 

The NOBLH iterative methods we use for the three- and six-variable models are ad hoc 
methods that arose from trying to glean more insight about our models’ behavior. Had 
we been interested solely in identifying an optimal solution, then an adaptive sequential 
procedure, such as response surface methodology [34], STRONG [35], or other simulation 
optimization approaches would require much less sampling. 

6.1.5 Robustness Analysis 

After determining the optimal variable values, our analysis turns to evaluating the robust¬ 
ness of the variables by evaluating power generation while incorporating deviations from 
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the optimal variable values. In the three-variable model, we find the separation between the 
frequency of excitation of the base (co) and natural frequency of the generator (co„) to be 
the most sensitive variable. While the designers of the piezoelectric generator should seek 
to minimize the generator’s modal damping ratio (£), the designer’s priority should be to 
precisely match the generator’s natural frequency with the frequency of excitation driving 
the generator. 

The precise matching of the base motion driving frequency {(Obase) to the modal short circuit 
net frequency (wsc) proves to be the most important factor in minimizing power loss during 
the robustness analysis of the six-variable model. Setting the electromagnetic coupling 
coefficient (A), load resistance (R), and net clamped capacity of piezoelectric material (Co) 
to their respective optimal value is a much lower priority for the piezoelectric generator 
designer. In fact, our robustness study concludes if the likelihood of negative variation is 
greater than or equal to positive variation for the design’s modal damping ratio, it is best 
to let £ vary from its optimal value while restraining the variation of other variables from 
their optimal values in order to produce the greatest mean power output. 

6.2 Further Studies 

Through the analysis of two mathematical models, we determine optimal values for each 
variable, given a practical variable range in order to maximize power generation of a piezo¬ 
electric generator. Additional studies are needed to propel our findings into practical appli¬ 
cation of piezoelectric materials for energy harvesting. 

6.2.1 PDE Solution 

To derive the nine variable model in Chapter 3, we use a one-mode expansion method. 
An alternative methods to derive models would be to use a partial differential equation 
expansion [36]. This mathematical derivation would provide an additional piezoelectrical 
generator model to conduct further optimization studies. 

6.2.2 Variable Range Adjustment 

We derive the variable ranges through research of materials specifications and previous 
experimental measurements. Personnel more experienced in vibration, piezoelectric mate¬ 
rials, and electrical engineering will have the ability to further adjust the ranges for each 
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variable. Using our methods of analysis and power generation models, it is possible to find 
updated optimal values for each variable. 

6.2.3 Simulation Optimization 

In future studies, an automated simulation optimization technique would be a better ap¬ 
proach for finding optimal nominal settings for the variables while keeping sampling re¬ 
quirements to a minimum. 

6.2.4 Robustness Analysis Expansion 

Our efforts study the robustness of the two power generation models based on both single 
and multiple variable manipulation. In both the three- and six-variable power manipulation, 
we assume the environment to be static. In other words, we assume the frequency of 
excitation (to) and the base motion driving frequency to be constant for the three- and 
six-variable models, respectively. This enables us to focus our robustness study on fewer 
variables. 

In practical applications, the frequency of the driving vibration may not prove constant. 
The vibration of an air conditioning duct may change based on temperature, humidity, fan 
speed, or a myriad other factors. Further robustness analysis focused on the base motion 
driving frequency may lead to a better piezoelectric generator power output in varying 
environmental conditions. 

6.2.5 Model Adjustment 

The two mathematical models we use do not serve as the only models to simulate piezoelec¬ 
tric power production. Further studies could include additional models or adjustments to 
the two models used in this study. The optimization and robustness analysis methodologies 
in this study could serve as a framework in the study of new models. 

6.2.6 Practical Experimentation 

We gear the efforts of this paper in order to find insights into maximizing piezoelectric 
generator power output. Using our findings, engineers can attempt to match our modeled 
results using a prototype piezoelectric generator designed with our optimal variable val¬ 
ues. Feedback from the results of the experimentation could serve to adjust our models or 
variable limits. 
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APPENDIX: MATLAB Codes 


A.l Three-Variable Orthogonality Computation Code 

%3 Variable Orthogonality Computation 
%Created April 2015 by Russell Nelson 

clear all 

%Import design points 
filename = 'Factorial3var.csv' ; 

%Rename file 
NOLH=csvread(filename) ; 

%Compute orthogonality for each variable 

avgomega = mean(NOLH(:, 1)) ; 

avgomega_n = mean(NOLH(:, 2)) ; 

avgzeta = mean(NOLH(: , 3)) ; 

numeratorl = 0; 

numerator2 = 0; 

numerator3 = 0; 

denomll = 0; 

denoml2 = 0; 

denom21 = 0; 

denom22 = 0; 

denom31 = 0; 

denom32 = 0; 

for i = 1:1000 

omegai = NOLH(i,l); 
omega_ni = NOLH(i,2); 

numeratorl = numeratorl + (omegai-avgomega)*(omega_ni-avgomega_n); 
denomll = denomll + (omegai-avgomega) A 2; 
denoml2 = denoml2 + (omega_ni-avgomega_n) A 2; 

end 

%Outputs correlation for omega to omega_n 
Ortho_omega_omega_n=numeratorl/sqrt(denomll*denoml2) 


83 





for i = 1:1000 

omega! = NOLH(i,l); 
zetai = NOLH (i,3); 

numerator! = numerator! + (omegai-avgomega)*(zetai-avgzeta); 
denom21 = denom21 + (omegai-avgomega)^2; 
denom22 = denom22 + (zetai-avgzeta)^2; 

end 

%Outputs correlation for omega to zeta 
Ortho_omega_zeta = numerator2/sqrt(denom21*denom22) 

for i = 1:1000 

omega_ni = NOLH(i,2); 
zetai = NOLH (i,3); 

numerator3 = numerator3 + (omega_ni-avgomega_n)*(zetai-avgzeta) 
denom31 = denom31 + (omega_ni-avgomega_n) * 2 ; 
denom32 = denom32 + (zetai-avgzeta)^2; 

end 

%Outputs correlation for omega_n to zeta 
Ortho_omega_n_zeta = numerator3/sqrt(denom31*denom32) 


A.2 Six-Variable Orthogonality Computation Code 

%6 Variable Orthogonality Computation 
%Created April 2015 by Russell Nelson 

clear all 

%Import design points 
filename = ' Factorial6var . csv' ; 

%Rename file 
NOLH=csvread(filename); 

%Compute orthogonality for each variable 
num=0; 
denl=0; 
den2=0; 
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for d = 2:6 

for i =1:15625 

Ai = NOLH (i,1) ; 

Bi = NOLH(i,d); 

avA = mean(NOLH(:,1)); 

avB = mean(NOLH(:,d)); 

nun = num + (Ai-avA)*(Bi-avB) 

deni = deni +(Ai-avA) A 2; 

den2 = den2 +(Bi-avB) A 2; 

end 

OrthoA(d)= num/sqrt(denl*den2); 

end 

num=0 ; 
denl=0; 
den2=0; 
for d = 3:6 

for i =1:15625 

Ai = NOLH (i,2) ; 

Bi = NOLH(i,d); 

avA = mean(NOLH ( :,2)); 

avB = mean(NOLH(:,d)); 

num = num + (Ai-avA)*(Bi-avB) 

deni = deni +(Ai-avA) A 2; 

den2 = den2 +(Bi-avB) A 2; 

end 

OrthoR(d)= num/sqrt(denl*den2); 

end 

num=0 ; 
denl=0; 
den2=0; 
for d = 4:6 

for i =1:15625 

Ai = NOLH(i,3); 

Bi = NOLH(i,d); 

avA = mean(NOLH(:,3)); 

avB = mean(NOLH(:,d)); 

num = num + (Ai-avA)*(Bi-avB) 
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deni = deni +(Ai-avA) A 2; 
den2 = den2 +(Bi-avB) A 2; 

end 

OrthoO(d)= num/sqrt(denl*den2) ; 

end 

num=0 ; 
denl=0; 
den2=0; 
for d = 5:6 

for i =1:15625 

Ai = NOLH (i,4) ; 

Bi = NOLH(i,d); 

avA = mean(NOLH ( :,4)); 

avB = mean(NOLH(:,d)); 

nun = num + (Ai-avA)*(Bi-avB); 

deni = deni +(Ai-avA) A 2; 

den2 = den2 +(Bi-avB) A 2; 

end 

OrthoC(d)= num/sqrt(denl*den2) ; 

end 

num=0; 
denl=0; 
den2=0; 
for d = 6:6 

for i =1:15625 

Ai = NOLH(i,5); 

Bi = NOLH(i,d); 

avA = mean(NOLH(:,5)); 

avB = mean(NOLH(:,d)); 

num = num + (Ai-avA)*(Bi-avB) ; 

deni = deni +(Ai-avA) A 2; 

den2 = den2 +(Bi-avB) A 2; 

end 

OrthoW(d)= num/sqrt(denl*den2) ; 

end 
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A.3 Euclidean Maximin Distance Code 


%Euclidean Maximin Distance Calculator 
%Created April 2015 by Russell Nelson 
clear all 

%Import factorial design points to use 
%as reference 

file = 'Factorial3var.csv' ; 
base = csvread(file); 

%Import Design Points 
filename = 'Factorial3var.csv' ; 

%Rename file 
A=csvread(filename) ; 

%Establish number of variables (n) 

%and number of design points (k) 

[n,k] = size (NOLH); 

d=0; 

i=l; 

%Normalizing vectors 
for p = 1:k 

avg = (max(base(:,p))-min(base(:,p)))/2+min(base(: , p)) 
for q = 1:n 

A(q,p) = (A(q,p)-avg)/avg; 

end 

end 

%Computing all distances 
for s = 1:n 

for t = 2:n 

if s == t 
else 

for j = 1:k 

d = d + abs(A(s,j)-A(t, j)) A 2 ; 
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end 

d = sqrt(d); 
distance (i) = d; 
d=0 ; 

i = i+1; 

end 

end 

end 

%Finding the Euclidean maximin distance 
min(distance) 


A.4 Three-Variable Model Optimization Code 
for Factorial Sampling 

%Optimal Power for Three Variable Model (Factorial Sampling Method) 
%Created January 2015 by Russell Nelson. 

%The following script evaluates the average power dissipated by the load 
%resistor of a piezoelectric generator. The equation, consisting of 3 
%independent variables, was derived by Professor Hong Zhou, 
clear all 

a=l; 
b=l; 
c=l ; 

%Set variable limits (range) 
omega_min_lim = 120*pi; 
omega_max_lim = 360*pi; 
omega_n_min_lim = 120*pi; 
omega_n_max_lim = 360*pi; 
zeta_min_lim = 0.005; 
zeta_max_lim = 0.02; 

%Set number of levels per variables 
levels = 100; 
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%Computes values of power for all combinations of 
%variables at all levels 

for omega = linspace(omega_min_lim,omega_max_lim,levels) 
b=l ; 

for omega_n = linspace(omega_n_min_lim,omega_n_max_lim,levels) 
c=l; 

for zeta = linspace(zeta_min_lim,zeta_max_lim,levels) 
P(a,b,c) = (zeta*((omega/omega_n) A 3)*(omega A 3))/ . . . 

((1-((omega A 2)/(omega_n A 2))) A 2+ ... 

(2*zeta*(omega/omega_n)) A 2) ; 
c=c+l; 

end 
b=b+l; 

end 
a=a+l; 

end 

%Vectorizes P matrix 
Pvec=P(:) ; 

%Find the maximum power in Pvec 
[v,I]=max(Pvec) 

%Determines the location of each variable level for maximum power 
[d, e,f]=ind2sub(size(P),I); 

%Determines the value for each variable at maximum power 
omega = linspace(omega_min_lim,omega_max_lim,levels); 
omegaopt=omega(d) 

omega_n = linspace(omega_n_min_lim,omega_n_max_lim,levels); 
omega_nopt=omega_n(e) 

zeta = linspace(zeta_min_lim,zeta_max_lim,levels); 
zetaopt=zeta(f) 

A.5 Three-Variable Model Optimization Code 
for NOBLH Sampling 


89 


%3 Variable Model Optimization Code for NOBLH Sampling 
%Created January 2015 by Russell Nelson 

%The following script evaluates the average power dissipated by the load 
%resistor of a piezoelectric generator. The equation, consisting of 3 
%independent variables, was derived by Professor Hong Zhou. 

clear all 

%Import design points 
filename = 'NOL3var.csv' ; 

%Rename file 
NOLH=csvread(filename); 

%Compute power for all design points 
for i =1:512 

omega=NOLH(i,1); 
omega_n=NOLH (i,2); 
zeta=NOLH (i,3); 

Power(l,i) = (zeta*((omega/omega_n) A 3)*(omega A 3))/ ... 

((1-((omega A 2)/(omega_n A 2))) A 2+(2*zeta*(omega/omega_n)) A 2); 

end 

%Identify maximum power and associated design point 
[x,y]=max(Power) 


A.6 Six-Variable Model Optimization Code 
for Factorial Sampling 

%Optimal Power for Six Variable Model (Factorial Sampling Method) 
%Created January 2015 by Russell Nelson 

%The following script evaluates the average power dissipated by the load 
%resistor of a piezoelectric generator. The equation, consisting of 6 
%independent variables, used was derived by Professor Hong Zhou. 
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clear all 


a=l 

b=l; 

c=l; 

d=l; 

e=l ; 

f=l; 

%Set variable limits(range) 

A_min_lim = 0.01; 

A_max_lim = 0.99; 

R_min_lim = .02; 

R_max_lim = 100; 

Obase_min_lim = 120*pi; 

Obase_max_lim = 360*pi; 

C0_min_lim = 0.00000001; 

C0_max_lim = 0.000001; 

WSC_min_lim = 120*pi; 

WSC_max_lim = 360*pi; 
zeta_min_lim = 0.005; 
zeta_max_lim = 0.02; 

%Set number of levels per variable 
levels=12; 

%Computes values of power for all combinations of 
%variables at all levels 

for A = linspace(A_min_lim,A_max_lim,levels) 
b=l; 

for R = linspace(R_min_lim,R_max_lim,levels) 
c=l ; 

for Obase = linspace(Obase_min_lim,Obase_max_lim,levels) 
d=l ; 

for CO = linspace(C0_min_lim,C0_max_lim,levels) 
e=l ; 

for WSC = linspace(WSC_min_lim,WSC_max_lim,levels) 
f=i; 

for zeta = linspace(zeta_min_lim,zeta_max_lim,levels) 
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P(a,b,c,d,e,f) =(A A 2*R*0base A 6*(1+CO A 2*0base A 2 *R A 2) A 2)/ ... 
(2*(((~l*Obase A 2+WSC A 2)*(l+CO A 2*Obase A 2*R A 2)+(C0*A A 2)*... 
Obase A 2*R A 2) A 2+(2*WSC*zeta*0base*(1+CO A 2*Obase A 2*R A 2)... 
+A A 2*Obase*R) A 2) ) ; 
f=f+1; 
end 
e=e+l; 
end 
d=d+l; 
end 
c=c+l; 
end 
b=b+l; 
end 
a=a+l; 
end 

%Vectorizes P matrix 
Pvec=P(:); 

%Finds the maximum power in Pvec 
[v,I] = max(Pvec) 

%Determines the location of each variable level for maximum power 
[j,k,l,m,n,o]=ind2sub(size(P) , I) ; 

%Determines the value for each variable at maximum power 
A = linspace(A_min_lim, A_max_lim,levels) ; 

Aopt = A(j) 

R = linspace(R_min_lim,R_max_lim,levels); 

Ropt = R(k) 

Obase = linspace(Obase_min_lim,Obase_max_lim, levels) ; 

Obaseopt = Obase(1) 

CO = linspace(C0_min_lim,C0_max_lim,levels); 

COopt = CO (m) 

WSC = linspace(WSC_min_lim,WSC_max_lim,levels) ; 

WSCopt = WSC(n) 

zeta = linspace(zeta_min_lim,zeta_max_lim,levels); 
zetaopt = zeta(o) 
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A.7 Six-Variable Model Optimization Code 
for NOBLH Sampling 

%6 Variable Model Optimization Code for NOBLH Sampling 
%Created January 2015 by Russell Nelson 

%The following script evaluates the average power dissipated by the load 
%resistor of a piezoelectric generator. The equation, consisting of 6 
%independent variables, used was derived by Professor Hong Zhou. 

clear all 

%Import design points 
filename = 'NOL6var.csv' ; 

%Rename file 
NOLH=csvread(filename); 

%Compute power for all design points 
for i =1:512 

A=NOLH(i,1); 

R=NOLH (i,2); 

Obase=NOLH (i,3); 

CO=NOLH(i,4); 

WSC=NOLH (i,5); 
zeta=NOLH(i,6); 

Power(1,i) =(A A 2*R*Obase A 6*(1+C0 A 2*Obase A 2*R A 2) A 2) / . . . 

(2*(((-l*Obase A 2+WSC A 2)*(l+CO A 2*Obase A 2*R A 2)+(C0*A A 2)*... 
Obase A 2*R A 2) A 2+(2*WSC*zeta*Obase*(1+C0 A 2*Obase A 2*R A 2)... 

+A A 2*Obase*R) A 2)) ; 

end 

%Identify maximum power and associated design point 
[x,y]=max(Power) 
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